HCl Accommodation, Dissociation, and Propensity for the Surface of

Article pubs.acs.org/JPCA

HCl Accommodation, Dissociation, and Propensity for the Surface of Water Collin D. Wick* Louisiana Tech University, P.O. Box 10348, Ruston, Louisiana 71270, United States ABSTRACT: A new reactive and polarizable molecular model was developed to describe HCl dissociation in liquid water and used to investigate HCl behavior at the air−water interface. It was found that the mechanism of HCl accommodation at the air−water interface began with its hydrogen pointing toward the water as it approached from the air. This was followed by dissociation into a contact ion pair once solvated at the air− water interface with the hydronium oriented more toward the air than the chloride on average. In comparison with NaCl, HCl showed some similar behavior in that its contact ion pair was stabilized at the air−water interface in comparison with the bulk. However, dissociated HCl had a greater propensity for the air−water interface than NaCl due to the fact that the hydronium ion was more surface active than sodium.

I. INTRODUCTION Interfacial environments hold unique attributes that distinguish them from bulk solvents. If certain species accumulate in the interfacial region, they may significantly alter chemical reactivity by catalyzing or inhibiting certain reactions. For instance, there has been recent debate on the concentration of acids and bases in the interfacial region.1−4 Many reactions are known to be catalyzed by the presence of an acid or a base, such as hydrolysis and transesterification.5,6 If it is determined that the interface can preferentially accommodate hydronium or hydroxide ions, this could potentially be exploited to increase the reactivity of these species. Additionally, atmospheric aerosols play an important role in the chemical transformations that shape atmospheric compositions.7−11 Any trace gas that comes in contact with an aerosol must first come in contact with its surface, making the interfacial region of enhanced importance. In particular, the aqueous concentration of HCl has been found to strongly influence the reactive uptake of trace gases.11 HCl itself is an important trace gas in the atmosphere that mainly forms from dechlorination of sea salt.12,13 HCl is known to react with aqueous aerosols to form many potential biproducts and is a potential source for rain acidity.14 HCl, as a strong acid, readily dissociates in water to hydronium and chloride ions. Less is known, though, about how this is affected by the presence of an interface and the dissociation mechanism itself. Small clusters with HCl have been extensively studied both experimentally and computationally.15−20 They bring insight into how more water molecules can lead to HCl dissociation, which includes contact ion pairs (CIPs) and solvent-separated ions pairs (SSIPs). More difficult to describe is how HCl behaves at an interface and comparing this behavior with bulk. QM/MM calculations have found that for HCl to dissociate a few water molecules have to cluster around it at the outer edge of a water surface.21 © 2013 American Chemical Society

What is less known, though, is how the interfacial region itself influences dissociation of a solvated HCl molecule. For nitric acid, it has been shown that hydrolysis is suppressed at the air− water interface in comparison with bulk.22,23 More controversial is the HCl propensity and orientation at the air−water interface. Addition of HCl to water decreases its surface tension,24 which, using the Gibbs surface tension equation, can be to related to a positive solute surface excess.25 In contrast, addition of NaCl to water increases its surface tension, which coincides with a negative solute surface excess.24 Taken together, this shows that at the same bulk concentration, the hydronium ion has a greater propensity for the air−water interface than sodium by around a factor of 1.53.24 Further support for this comes from how the addition of HCl changes the mobility of aqueous droplets in capillary electrophoresis measurements, which suggests that the surface of water droplets becomes more positively charged with addition of HCl.1 In contrast, addition of HCl increases the surface potential, which suggests that negative charge accumulates toward the surface of water with addition of HCl.26 Computational studies of the air−water interface disagree with this for mixtures of hydronium and chloride,27−29 with debate over the degree of hydronium propensity for the air−water interface. Interestingly, though, molecular HCl at the air−water interface has been found to orient with its hydrogen pointing toward the bulk, in contrast to what was found for the dissociated state.21 If hydronium ions preferentially reside at the air−water interface, as suggested, it may significantly enhance the reactivity of different species present there. For instance, in the reactive uptake of trace gases, they would come in contact with a more Received: August 22, 2013 Revised: October 25, 2013 Published: October 29, 2013 12459

dx.doi.org/10.1021/jp4084212 | J. Phys. Chem. A 2013, 117, 12459−12467

The Journal of Physical Chemistry A

Article

Figure 1. Schematic of the implementation of the MS-EVB model for a HCl cluster. NS

acidic environment than is present in bulk, potentially significantly affecting their uptake mechanism. Classical molecular simulations of HCl solutions treat hydronium and chloride ions as completely separate molecules.28,29 However, it is well known that hydronium ions can share their charge with other water molecules.30 Furthermore, to have a description of HCl dissociation, a reactive molecular model is required. Ab initio methods are often preferred to investigate reactive systems because they do not require significant molecular model development. HCl aqueous systems have been investigated with ab initio molecular dynamics (AIMD),31 along with QM/MM calculations.21 AIMD methods with some of the most powerful supercomputers can, in principle, be used to investigate large aqueous systems with interfaces.32 However, a computationally less expensive methodology is sought. The multistate empirical state valence-bond (MS-EVB) method is a good balance between computational power and the ability to describe the necessary physics for proton solvation and transfer.33−37 The MS-EVB model treats a proton as being in multiple states, scales linear with the number of states treated, and is much less computationally expensive than AIMD calculations. This method has been used to bring insight into the behavior of hydronium ions in a variety of environments,33−38 including its role in the dissociation of a weak acid.36 As a consequence, it is an attractive method for investigation of HCl dissociation in that it both is reasonably efficient and can describe acid dissociation. This work describes a computational investigation of HCl dissociation at the air−water interface using the MS-EVB model with inclusion of polarizable interactions. Comparisons with NaCl dissociation are also made to discern the differences between the interfacial behavior of acid and salt solutions. A new HCl model was developed based off previously published hydronium and chloride ion models.39,40

HEVB =

NS

∑ ∑ ⟨i|Ĥ |j⟩ i=1 j=1

(1)

in which NS encompasses all potential hydronium states. The potential hydronium states are determined as follows. A hydronium donates up to three hydrogen bonds, and each of the water molecules that receive one is considered to have hydronium character. These first solvation shell water molecules can donate up to two more hydrogen bonds a piece to other water molecules, and each of these second solvation shell water molecules are considered to have hydronium character as well. In eq 1, Ĥ is the Hamiltonian for individual states, each with the hydronium ion present on a different oxygen. Classical hydronium and water interaction potentials are used to calculate diagonal terms in a NS × NS matrix with the hydronium ion present on a different oxygen in each state with the remaining water oxygens treated as part of a water molecule for each state. Nondiagonal terms in this matrix are parametrized to reproduce the potential energy surfaces (PESs) of different ab initio clusters, ranging from H5O2+ to H9O4+ in size. The energy of the system is calculated by finding the lowest energy eigenvalue of the matrix.37 The location of the hydronium ion is determined to be the state with the highest squared eigenvector and is used to determine the hydronium position in each step. For the investigation of strong acids and bases, it has generally been assumed that the dissociated state is enough to describe the systems due to the rarity of the observance of molecular species.27,42 For HCl, its pKa is very negative, around −7,43 making it expected that all solvated species can be treated as dissociated. However, to describe the sorption of gaseous HCl onto the surface of water, which is relevant for atmospheric uptake, molecular HCl needs to be taken into account. Figure 1 shows a schematic of the implementation of the MS-EVB model for HCl in a cluster of four molecules. In this example, the state with the highest eigenvector has the hydronium ion in the middle (state 1). If only water and hydronium states were described, states 1, 4, 5, and 6 would be all that is represented. However, to account for the potential bound HCl, states 2 and 3 need to be included also. It should be noticed that the chloride can potentially bind to more than

II. MOLECULAR MODEL DETAILS The hydronium and water models are detailed elsewhere.39 The model utilizes the MS-EVB framework, which has been employed for many different systems,33−35,37,38,41 including strong acids and bases.27,42 Briefly, the MS-EVB model uses a Hamiltonian to describe each potential hydronium state 12460



Article

literature.50 Also, α = 1.913 and r0 = 1.27332 Å, which were parametrized to reproduce the HCl stretch energy with respect to distance from MP2 calculations with the aug-cc-pvtz basis set. U0 was set to −165.25 kcal/mol to give the same difference between [HCl + H2O separated] and [H3O+ + Cl− separated] clusters in comparison with ab initio calculations (MP2 with aug-cc-pvtz basis set). Sodium and chloride potentials had charges and point polarizabilities located on their atomic positions, taken from gas-phase ab initio calculations.51,52 They have integer charges, and their exp-6 σ and ε parameters are given in Table 1. These were parametrized to reproduce the free energy of solvation calculated via the staged free energy perturbation (FEP) method.53 Further details of their parametrization are given in other work.54 II.B. Nondiagonal Terms. In general, nondiagonal terms are present to link states that share a hydrogen. For instance, in Figure 1, states 1 and 4 share a hydrogen (the one just below the center oxygen), states 4 and states 5 share one, and states 1 and states 6 do too. There are cross terms to address all of these (h14, h45, and h16). These were parametrized to reproduce MP2 calculations of H3O+ with up to three water molecules present and were functions of rOO and rOH between hydrogen binding species, where the H was bisecting the two Os.39 The form of this potential was different depending on the number of water molecules the hydronium ion donated hydrogen bonds to. Thus, different potentials were used if a hydronium ion donated one hydrogen bond, two, and three with switching functions present to make sure transitions from one potential to another were smooth

one hydrogen and is why there are two states shown in Figure 1 describing molecular HCl. II.A. Diagonal Terms. As with the situation of hydronium in water, each of these states is calculated with a classical molecular model. The hydronium and water models are described elsewhere,39 while a new model was developed for HCl. The repulsive and van der Waals interactions were described with a single Buckingham exponential-6 (exp-6) potential site located on the heavy atom positions 6 ⎡6 ⎛ ⎡ r ⎤⎞ ⎛ σ ⎞ ⎤ Uvdw = ε⎢ exp⎜λ⎢1 − ⎥⎟ − ⎜ ⎟ ⎥ σ ⎦⎠ ⎝ r ⎠ ⎦ ⎣λ ⎝ ⎣

(2)

with sigma and epsilon values given in Table 1 and λ = 13.5. Unlike interactions are treated with the standard Lorentz− Table 1. Parameters Used To Describe the Classical Interactions of HCl and Cl− atom

σ (Å)

ε (Å)

q (e)

α (Å3)

Cl− (H)−Cl H−(Cl)

4.89 0 3.96

0.1 0 0.6367

−1.0 0.23 −0.23

5.482 0.0 1.91

Berthelot combining rules (arithmetic mean for σ and geometric mean for ε). For HCl, the values were estimated by taking Lennard−Jones parameters from a previously developed HCl model44 and scaling them in a similar fashion as for the water oxygen going from a Lennard−Jones to an exponential-6.39,41 It should be noted that most of the interactions between HCl and water are described by a combination of diagonal terms and nondiagonal terms (described in the next section) in the MS-EVB matrix, which are parametrized to ab initio calculations. The HCl chloride atom had a point polarizability on its atomic position with a value taken from the work of Applequist.45 Partial charges were placed on the hydrogen and chloride atomic positions given in Table 1. These were chosen based on ESP charges with 0.1 Å grid spacing calculated via a MP2 calculation of HCl with the aug-cc-pvtz basis. Short-ranged electrostatic interactions were damped via the method described by Thole46 and further described by other workers.47,48 Specific implementation is given in detail elsewhere.39 A Morse potential was used to describe HCl bond stretching of the form UBond(r ) = D[1 − e−α(r − r0)]2 + U0

Urswitch =

1 {1 − tanh[40(rXH − r0)]} 2

(5)

where rXH is the distance between oxygen or chloride and hydrogen with r0 = 2.1 Å for oxygen−hydrogen distances and a longer r0 = 2.6 Å for HCl distances. These are designed so the cross terms smoothly go to zero when an oxygen−hydrogen distance hits 2.3 Å, which is the cutoff for hydronium−water interactions to be considered for the MS-EVB potential, and 2.8 Å for HCl distances. In Figure 1, states 2 and 3 are used to describe molecular HCl, which required additional cross terms (h12 and h16) to describe. These are simply a function of rHCl, which was found to be adequate to describe HCl dissociation in water. Figure 2 gives the three clusters used to parametrize the HCl nondiagonal terms. For the HCl cluster with one water molecule present (cluster 1), it can be observed that HCl is still strongly associated while cluster 4a shows a CIP and 4b shows a SSIP, both with four water molecules present in addition to

(3)

where D = 106.2 kcal/mol, which was extracted from the experimental D0 value49 plus 1/2hν0, with ν0 taken from the

Figure 2. Snapshots of structures used to calculate the energies for the 1, 4a, and 4b clusters. 12461



Article

HCl. The potential energy as a function of HCl distance was calculated from MP2 calculations with the aug-cc-pvtz basis set for cluster 1 and the aug-cc-pvdz for the clusters 4 and 4a with the NWChem computational package using counterpoise corrections.55,56 For cluster 1, specifically, it is described with two states, one with HCl and H2O (state 1) and another with Cl− and H3O+ (state 2). Each of the states are described via their classical models to account for the diagonal terms (h11 and h22), while h12 was parametrized to reproduce the potential energy as a function of rHCl. Instead of using a specific functional form, a cubic spline interpolation was used to parametrize h12 as a function of rHCl. The spline is based on a set of points for h12 spaced in 0.1 Å increments for rHCl, and the values are given in Table 3. For rHCl = 2.1 Å and greater, h12 = 0. Because of this, the switching function described in eq 5 is not necessary for cluster 1 but is used for larger clusters, as the water−hydronium cross terms are affected by the number of hydrogen bonds a hydronium donates (as discussed previously). Table 2 gives the ab initio and MS-EVB model energies for different clusters examined for this work. In addition to the

Figure 3. Potential energy as a function of HCl distance for the 1 (bottom) and 4a (top) clusters. Dotted lines represent the results for inclusion of NQE with the black lines/symbols representing the PIMC results and the red line/symbols the model fit to PIMC.

Nuclear quantum effects (NQE) may play an important role in HCl dissociation because of the importance of the very light H atom in this reaction. There are many ways to handle NQE. One of them is to run centroid molecular dynamics,57 which increase the computational cost by an order of magnitude. Another method is to create an effective potential based on the centroid of the results of path integral Monte Carlo (PIMC). This has been done multiple times before for hydronium− water and hydroxide−water interactions.38,39,41,58 In the current work, PIMC simulations were carried out for the 1 cluster and a new potential with respect to rHCl was fit to the distance between the centroid of the hydrogen and chloride atoms. As described previously, a cubic spline interpolation was used to describe h12 as a function of HCl distance based on a parametrization in 0.1 Å increments. These values are given in Table 3: hij (NQE). The new potential energy extracted from PIMC simulations of cluster 1 with respect to rHCl is given in Figure 3, along with results for our model after refitting to this, showing little change.

Table 2. Energy Comparisons Between the ab Initio and the MS-EVB Model for Different Gas-Phase Clusters Eab initio (kcal/mol)

Emodel (kcal/mol)

0.0 −165.25 −171.2 −199.42 −200.42

0.0 −165.25 −171.6 −199.56 −200.16

Cl− + H3O+ (separated) HCl + H2O (separated) HCl + H2O (bound, 1) HCl + 4H2O (close, 4a) HCl + 4H2O (far, 4b)

Table 3. Values Used for the Cubic Spline Interpolation for the Cross Terms for HCl Binding with Water rHCl

hij (ab initio)

hij (NQE)

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0

18.0 18.0 18.0 16.0 14.0 14.5 10.0 7.0 1.0 1.0

0.0 14.0 16.0 6.0 6.0 9.4 7.6 6.3 4.4 0.0

III. SIMULATION DETAILS All simulations included 1000 water molecules and were carried out at 298 K utilizing the Berendsen thermostat to control the temperature.59 We carried out simulations for two types of systems: bulk and interfacial. Bulk systems included a box with dimensions of approximately 28 Å × 28 Å × 38 Å with its pressure set to 1 atm controlled with the Berendsen barostat. For interfacial systems, a constant volume simulation was carried out in which the z dimension was elongated to 80 Å. Due to this elongation, water occupied approximately one-half the simulation box and air occupied the other half. As a consequence, two air−water interfaces formed. The Gibbs dividing surface (GDS) was determined by fitting a tanh function to the atomic density of water with the GDS representing the position of the half height of the function.60 A potential truncation of 9 Å was enforced for van der Waals interactions with analytical tail corrections employed, and the particle mesh Ewald summation was used for electrostatic interactions.61 The potential of mean force (PMF) with respect to the HCl interatomic position for both positions and the chloride position was mapped out for the systems. This was carried out utilizing umbrella sampling,62 which uses a pseudopotential

clusters shown in Figure 2, the energies of isolated molecules are shown. The zero of energy was chosen to be the separated Cl− and H3O+ ions. The other cluster with separated HCl and H2O did not include any contribution from the off-diagonal MS-EVB terms but simply were a result of the classical models. The three clusters shown in Figure 2, with energies given in Table 2, show very good agreement between simulations and ab initio calculations. Figure 3 shows a comparison between our model and ab initio calculations for the potential energy as a function of HCl distance for clusters 1 and 4a. The agreement between ab initio and our model is quite good for the 4a cluster while not so good for the 1 cluster but still reasonable. The reason for putting greater emphasis on reproducing the energy of the 4a cluster is that nearly all simulations will include HCl in contact with many water molecules. 12462


The Journal of Physical Chemistry A u(z) = kr(rXCl − r0)2 + kz(zA − z 0)2

Article

(6)

where rXCl represents the HCl or NaCl distance and zA is the distance from the water slab center of mass minus the average position of the GDS, making a value of 0 represent the position at the average GDS. For NaCl systems, A is the sodium ion position, and for HCl systems, it represents the hydronium ion position with the highest MS-EVB eigenvector. The force constants, kr and kz, are set to 0 or 2.0 kcal mol−1 Å2 depending on the type of PMF being carried out. The PMFs are with respect to either rXCl or zA for this work. When creating a PMF with respect to one coordinate, such as rXCl, multiple simulations are carried out, each with different values of r0. The umbrella potential will keep the coordinate centered around r0 for that specific simulation. Then, the PMFs extracted from each simulation are matched in overlapping regions to construct a PMF over a wide range of rXCl values for this example. For all PMF calculations with respect to rHCl, 200 ps of simulation time were carried out with 10 different r0 values used (0.5, 1.0, 1.5, 1.7, 2.0, 3.0, 4.0, 5.0, 6.0, and 7.0 Å). For PMF calculations with respect to rNaCl, 400 ps of simulation time for each r0 values were used (2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0 Å). All PMF calculations with respect to zA included values of z0 equal to −8, −6, −4, −2, −1, 0, 1, and 2 Å with 0 representing the GDS and negative values represent the liquid phase. All uncertainties in the PMFs were calculated by splitting the simulation results into four separate blocks and calculating the standard error of the mean.

Figure 5. PMF for HCl dissociation with respect to HCl interatomic distance in bulk water and at the air−water interface.

by the location of HCl. There are noticeable local minimums at 1.3 Å, which is clearly associated with molecular HCl. After that there appears to be a large peak, followed by a minimum around 2.0 Å, representing the CIP. The free energy minimum for the CIP is shifted to shorter distances for the interfacial systems, showing a stabilization of the CIP at the air−water interface. This is similar to what was found for NaCl. As with NaCl, chloride polarizability induces a higher chloride induced dipole at the air−water interface than in the bulk.65 These higher induced dipoles should promote stronger ion pairing at the interface and are likely the reason that the CIP is stabilized at the interface for both NaCl and HCl. Interestingly, the differences in the minimums between molecular HCl and the CIP are similar at both the bulk and at the interface. One may think molecular HCl would be stabilized more at the interface, but this is offset by the stabilization of the CIP there also. IV.B. HCl PMF with Respect to zA. Figure 6 gives the PMF with respect to the z position for hydronium across the air−

IV. RESULTS AND DISCUSSION IV.A. PMF With Respect to rXCl. The PMF with respect to rNaCl is given in Figure 4 in the bulk and at the air−water

Figure 4. PMF for NaCl dissociation with respect to NaCl interatomic distance in bulk water and at the air−water interface.

Figure 6. PMF across the air−water interface for the chloride ion for HCl pairs or the hydronium cation center of mass (see text for more info). Zero represents the GDS.

interface. For calculations at the air−water interface, z0 was set to 19 Å from the water liquid center of mass, approximately equal to the GDS (which was 19.2 Å from the water center of mass) with kz = 2.0 kcal/mol. NaCl is more strongly associated at the air−water interface than in the bulk with a noticeably lower free energy minimum at rNaCl = 2.8 Å. This is consistent with previous work that found the same for another polarizable NaCl model.63,64 Figure 5 also gives the free energy as a function of HCl interatomic distance with z0 set at 19 Å from the water liquid center of mass for the interfacial system. The H in HCl position was assigned based on the hydrogen of the hydronium ion with the highest eigenvector (state 1 in Figure 1) that was closest to the chloride position. The PMF’s overall shape and structure does not appear to be significantly affected

water interface for different HCl separations. It should be noted that the distance shown is the r0 values in eq 6, so the actual distance was not fixed but was a distribution close to that value. These distances were chosen to relate to the minimums in the PMF with respect to rHCl, corresponding with molecular HCl (r = 1.2 Å), the CIP (r = 2.2 Å), and SSIP (r = 4.2 Å). The z position is with respect to the average hydronium position which was calculated by determining its expectation value for a given configuration 12463


The Journal of Physical Chemistry A NS

x EVB =

Article

NS

potential as a function of z position is calculated from partial charges and induced dipoles69−71

∑ ∑ ⟨i|x|j⟩ i=1 j=1

(7)

Δϕq(z) = ϕq(z) − ϕq(z 0)

which results in the position being somewhat delocalized. The only exception is for molecular HCl, in which the chloride position is shown due to the fact that hydronium states are not populated to a significant degree. Figure 6 shows that molecular HCl has a significant free energy minimum at the interface, corresponding to the region where it would be partially solvated, followed by a free energy barrier, and then the PMF levels out as HCl moves further toward the bulk. Many trace gases, such as the hydroxyl radical and ozone, have a free energy minimum at the air−water interface.66,67 Unlike these described trace gases, though, HCl readily dissociates into hydronium and chloride ions. This requires an understanding of the CIP and more dissociated states to describe HCl solvation. The CIP for HCl shows no significant free energy minimum at the air− water interface, becoming lower as the hydronium ion approaches water until it levels off around 1 Å from the GDS toward the bulk. However, it can be observed that the hydronium ion has a rather low free energy in regions of positive z position, i.e., with its position more toward the air than the GDS. This is consistent with the picture of the hydronium ion present on the outer surface of water.27 For HCl in a SSIP, it has a small free energy minimum near the GDS, only slightly below the CIP. This slightly more prominent minimum for the SSIP is likely due to the hydronium interaction with chloride for the CIP, which will be discussed later. IV.C. HCl Orientation Near the Interface. Figure 7 gives both the hydronium and the chloride PMFs across the air−

=− Ez(z) =

1 ε0

∫z

z 0

∫z

⎛ ⎞ 1 ⎜E(z′) + ⟨ρμind (z′)⟩⎟dz′ ε0 ⎝ ⎠

(8)

z

⟨ρq (z′)⟩dz′ 0

(9)

where ρq is the charge density as a function of z position and ρμ is the dipole density as a function of position. There are two potential ways to evaluate the charge and dipole densities. One is to use the state which has the highest eigenvector, giving the hydronium and chloride ions single positions per configuration. This will be denoted the “classical” position. However, the actual charge and dipole on a particular atom can be calculated by its expectation value NS

AEVB =

NS

∑ ∑ ⟨i|A|j⟩ i=1 j=1

(10)

where A can either be the charge of an atom or its induced dipole for a particular configuration. The hydronium and chloride charges are delocalized to a degree, which may influence the electrostatic potential. Both the real and the classical electrostatic potentials were calculated from the PMF calculations with respect to z in which z0 = −4 Å from the GDS. Figure 8 gives the electrostatic potentials calculated as a

Figure 8. Electrostatic potential across the air−water interface from the CIP and SSIP PMF calculations with z0 = −4 Å, along with that for pure water. Figure 7. PMF across the air−water interface for both hydronium and chloride ions for HCl, CIP, and SSIP.

function of z position across the air−water interface, along with that for pure water calculated in previous work.72 The classical electrostatic potential becomes more negative with addition of HCl, which is consistent with the observation that the hydronium has a greater propensity for the air−water interface than the chloride ion, but in contrast to experimental electrostatic potential measurements.26 On the other hand, the regular electrostatic potential is more positive than that of pure water, showing agreement with experiment. Previous work has found that the hydronium electrostatic charge delocalizes toward the water bulk near the interface,27,39,73 which should increase the electrostatic potential across the interface. This may be a consequence of image charge repulsion with the interface74,75 or could be because of the electrostatic potential of pure water.76 The shoulder for the regular electrostatic potential for the CIP and SSIP at −3 Å that is not present for the classical potentials is most likely due to a shift of positive

water interface for the CIP and SSIP hydronium−chloride pairs. It is clear from the PMFs that hydronium is more likely to be found toward the air in comparison with chloride. The hydronium ion has its free energy increasing at 1 Å from the GDS toward the air and chloride’s increases 1 Å from the GDS toward the liquid bulk. This is consistent with the zeta potential changes with addition of HCl to water, which suggests that a double layer forms with the surface being positively charged.1,68 However, this is in contrast to the interpretation of the change in surface potential with addition of HCl,26 consistent with previous simulations of HCl.29 To better quantify the perceived discrepancy between simulation and experiment for the surface potential, the surface potential was calculated for this work. The electrostatic 12464



Article

IV.E. Comparison of NaCl and HCl Propensity for the Air−Water Interface. To better compare the differences between hydronium and sodium, Figure 10 gives the PMFs for

charge toward the bulk and a shift of negative charge toward the interface, creating a region with little to no net dipole with respect to the normal of the interface. This also shows the importance of including charge delocalization in the description of acids, as they can qualitatively influence interfacial properties. It has been suggested though that one should be careful in interpreting electrostatic potentials across the air−water interface, especially for pure water,77 but the change in potential with addition of inorganic salts71 and now addition of HCl shows that qualitative agreement with experiment is possible for this. Furthermore, one should caution that the results are from a single point used to construct the PMFs and direct comparisons with experiment require concentrated solutions without any umbrella sampling. IV.D. Mechanism for HCl Uptake. For the uptake of HCl into water, it starts out as molecular HCl, so the orientation of molecular HCl while near the interfacial region is important to describe this. To estimate the orientation molecular HCl as it comes in contact with the water surface, its average orientation was calculated from configurations used to sample z0 = 0.0, 1.0, and 2.0 Å (i.e., those located more toward the air than the GDS of water). From these, the average HCl vector with the z axis was found to have cos θ = −0.22 ± 0.14, −0.38 ± 0.11, and −0.76 ± 0.09 for z0 = 0.0, 1.0, and 2.0 Å, respectively, with negative values corresponding to the HCl hydrogen pointing toward the water bulk. This shows that as molecular HCl comes in contact with the outer edge of the water surface its hydrogen points toward the water surface and as it reaches the interface, it aligns more parallel to the interface but still with a preference to have its hydrogen pointing toward the bulk. Using the computational results, the mechanism for uptake of HCl onto water can be developed with Figure 9, showing

Figure 10. PMF across the air−water interface for the chloride ion in NaCl and HCl.

sodium from NaCl and the hydronium from dissociated HCl across the air−water interface. The CIP of NaCl has a free energy minimum for sodium positions between 2 and 4 Å from the GDS toward the bulk but becomes positive when it approaches the GDS. This compares with the CIP of HCl, which shows no free energy minimum for hydronium, but its free energy does not become positive until it is past the GDS toward the air. In general, this points to hydronium being more surface active than sodium, which is consistent with experimental surface tension measurements.24 In contrast, there appears to be a greater stabilization of sodium for a specific position near the interface that is not shared by hydronium. The SSIP of NaCl is similar to that of its CIP, but there is no free energy minimum near the air−water interface or it is much smaller than for the CIP. This is also in contrast to HCl, which shows a greater free energy minimum for the SSIP than for the CIP. The following is apparent. NaCl being in a CIP stabilizes it at the air−water interface to a greater degree than the CIP for HCl. For NaCl, chloride’s higher induced dipole when it is in a CIP is a likely reason it is stabilized at the air−water interface.64 HCl though has the hydronium oriented toward the air not the chloride, and inducing a dipole on the hydronium would cause its oxygen to behave as it was more negatively charged in the direction opposite of its hydrogen atoms. This would most likely enhance interactions with water hydrogens, reducing its propensity for the air−water interface in comparison with the SSIP. To better quantify the predicted differences between sodium and hydronium interfacial concentrations, the PMFs were converted into partition coefficients using the relation K(z) = exp(−W(z)/RT). Next, the partition coefficients were integrated as a function of z position, subtracting 1 only for z < 0. The integrated partition coefficient for the HCl hydronium was then divided by that for the NaCl sodium to estimate the ratio in surface excesses. The ratio arrived at for this was 1.43 with an estimated uncertainty of 0.2, which compares to the experimental value of 1.53, giving good agreement. It is clear that the hydronium has a greater propensity for the air−water interface than sodium but of interest is that the mechanism for hydronium interfacial behavior is much different than sodium. Sodium clearly has its position shifted more toward the liquid bulk, and hydronium clearly presides more toward the air. This is likely because of the fact that hydronium can delocalize its charge and that its oxygen does not strongly interact with water

Figure 9. Snapshots of configurations of HCl showing a proposed mechanism for HCl accommodation on the air−water interface.

snapshots illustrating this. First, molecular HCl comes in contact with the surface of water with its hydrogen pointing toward the liquid (Figure 9A). Next, as HCl approaches the GDS, its orientation shifts to be more parallel to the surface. There is a free energy minimum in this region for molecular HCl, followed by a barrier for further penetration into water. As a consequence, this is by far the most likely region for HCl dissociation to take place, which is consistent with previous work.21 Next, dissociation of HCl occurs to form a CIP, which is fairly stable in the interfacial region. Moreover, the CIP exists with the hydronium ion oriented toward the air and the chloride ion oriented toward the bulk (Figure 9B). However, there is not a free energy minimum associated with this state with respect to the z coordinate, which shows that once HCl dissociates it readily moves into the water bulk as either a CIP or a SSIP (Figure 9C). 12465



Article

(13) Graedel, T. E.; Keene, W. C. Global Biogeochem. Cycles 1995, 9, 47−77. (14) Sanhueza, E. Tellus, Ser. B: Chem. Phys. Meteorol. 2001, 53, 122− 132. (15) Arillo Flores, O. I.; Bernal-Uruchurtu, M. I. J. Phys. Chem. A 2010, 114, 8975−8983. (16) Ando, K.; Hynes, J. T. J. Phys. Chem. B 1997, 101, 10464− 10478. (17) Milet, A.; Struniewicz, C.; Moszynski, R.; Wormer, P. E. S. J. Chem. Phys. 2001, 115, 349−356. (18) Weimann, M.; Fárník, M.; Suhm, M. A. Phys. Chem. Chem. Phys. 2002, 4, 3933−3937. (19) Odde, S.; Mhin, B. J.; Lee, S.; Lee, H. M.; Kim, K. S. J. Chem. Phys. 2004, 120, 9524−9535. (20) Gutberlet, A.; Schwaab, G.; Birer, Ö .; Masia, M.; Kaczmarek, A.; Forbert, H.; Havenith, M.; Marx, D. Science 2009, 324, 1545−1548. (21) Ardura, D.; Donaldson, D. J. Phys. Chem. Chem. Phys. 2009, 11, 857−863. (22) Wang, S.; Bianco, R.; Hynes, J. T. J. Phys. Chem. A 2009, 113, 1295−1307. (23) Shamay, E. S.; Buch, V.; Parrinello, M.; Richmond, G. L. J. Am. Chem. Soc. 2007, 129, 12910−12911. (24) Pegram, L. M.; Record, M. T., Jr. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 14278−14281. (25) Atkins, P.; de Paula, J. Atkins’ Physical Chemistry, 5th ed.; University Press: Oxford, U.K., 1994. (26) Randles, J. E. B. Phys. Chem. Liq. 1977, 7, 107−179. (27) Iuchi, S.; Chen, H.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2009, 113, 4017−4030. (28) 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. (29) Ishiyama, T.; Morita, A. J. Phys. Chem. A 2007, 111, 9277−9285. (30) Eigen, M. Angew. Chem., Int. Ed. Engl. 1964, 3, 1−19. (31) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. Nature 1999, 397, 601−604. (32) Mundy, C. J.; Kuo, I. F. W.; Tuckerman, M. E.; Lee, H. S.; Tobias, D. J. Chem. Phys. Lett. 2009, 481, 2−8. (33) Swanson, J. M. J.; Maupin, C. M.; Chen, H.; Petersen, M. K.; Xu, J.; Wu, Y.; Voth, G. A. J. Phys. Chem. B 2007, 111, 4300−4314. (34) Voth, G. A. Acc. Chem. Res. 2006, 39, 143−150. (35) Day, T. J. F.; Soudackov, A. V.; Cuma, M.; Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 2002, 117, 5839−5849. (36) Cuma, M.; Schmitt, U. W.; Voth, G. A. J. Phys. Chem. A 2001, 105, 2814−2823. (37) Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 1999, 111, 9361− 9381. (38) Brancato, G.; Tuckerman, M. E. J. Chem. Phys. 2005, 122, 224507. (39) Wick, C. D. J. Phys. Chem. C 2012, 116, 4026−4038. (40) Cummings, O. T.; Wick, C. D. J. Chem. Phys. 2013, 139, 064708−064709. (41) Wick, C. D.; Dang, L. X. J. Phys. Chem. A 2009, 113, 6356− 6364. (42) Wick, C. D.; Dang, L. X. J. Chem. Phys. 2010, 133, 024705− 024708. (43) Wade, L. G. Organic Chemistry, 6th ed.; Prentice-Hall: Upper Saddle River, NJ, 2005. (44) Huang, Y. L.; Heilig, M.; Hasse, H.; Vrabec, J. AlChE J. 2011, 57, 1043−1060. (45) Applequist, J.; Carl, J. R.; Fung, K. K. J. Am. Chem. Soc. 1972, 94, 2952−2960. (46) Thole, B. T. Chem. Phys. 1981, 59, 341−350. (47) Burnham, C. J.; Li, J.; Xantheas, S. S.; Leslie, M. J. Chem. Phys. 1999, 110, 4566. (48) Sala, J.; Gurdia, E.; Masia, M. J. Chem. Phys. 2010, 133. (49) Martin, J. D. D.; Hepburn, J. W. J. Chem. Phys. 1998, 109, 8139−8142. (50) McQuarrie, D. A. Statistical Mechanics; Harper Collins: New York, 1973.

molecules. In other words, hydronium can behave in an amphiphilic manner, which has been promoted in previous work.33,34

V. CONCLUSIONS A new dissociable HCl molecular model was developed to understand its behavior at the air−water interface. The model was parametrized to ab initio calculations and found to give a dissociation constant in reasonable agreement with experiment. The HCl contact ion pair was found to be somewhat stabilized at the air−water interface in comparison with bulk, similar to what was found for NaCl. Molecular HCl had a very strong propensity for the air−water interface, while the contact ion pair had no free energy minimum there. The more dissociated solvent-separated ion pair had a small free energy minimum near the air−water interface. Comparing HCl with NaCl found that the hydronium had around 1.4 times the propensity for the air−water interface than sodium, consistent with experimental surface excess estimates. Finally, the mechanism for HCl accommodation at the air−water interface begins with molecular HCl having its hydrogen pointing toward the water surface until it comes in contact with water, in which its orientation becomes parallel to the interface. HCl then dissociates to a contact ion pair, followed by further dissociating or movement toward the water bulk, in which the hydronium is preferentially oriented toward the air in comparison with the chloride ion. NaCl ions have the opposite orientational preference.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was funded by a grant from the National Science Foundation (NSF EPS-1003897). Calculations were carried out using the resources from the Louisiana Optical Network Initiative (LONI).

■

REFERENCES

(1) Creux, P.; Lachaise, J.; Graciaa, A.; Beattie, J. K.; Djerdjev, A. M. J. Phys. Chem. B 2009, 113, 14146−14150. (2) Beattie, J. K.; Djerdjev, A. M.; Warr, G. G. Faraday Discuss 2008, 141, 31−39. (3) Buch, V.; Milet, A.; Vacha, R.; Jungwirth, P.; Devlin, J. P. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 7342−7347. (4) Petersen, P. B.; Saykally, R. J. Chem. Phys. Lett. 2008, 458, 255. (5) Polydore, C.; Roundhill, D. M.; Liu, H. Q. J. Mol. Catal. A: Chem. 2002, 186, 65−68. (6) Yu, Y.; Lou, X.; Wu, H. Energy Fuels 2008, 22, 46−60. (7) Prather, K. A.; et al. Proc. Natl. Acad. Sci. 2013, 110, 7550−7555. (8) Nissenson, P.; Knox, C. J. H.; Finlayson-Pitts, B. J.; Phillips, L. F.; Dabdub, D. Phys. Chem. Chem. Phys. 2006, 8, 4700−4710. (9) Knipping, E. M.; Lakin, M. J.; Foster, K. L.; Jungwirth, P.; Tobias, D. J.; Gerber, R. B.; Dabdub, D.; Finlayson-Pitts, B. J. Science 2000, 288, 301−306. (10) Swartz, E.; Shi, Q.; Davidovits, P.; Jayne, J. T.; Worsnop, D. R.; Kolb, C. E. J. Aerosol Sci 1998, 29, S985. (11) Davidovits, P.; Kolb, C. E.; Williams, L. R.; Jayne, J. T.; Worsnop, D. R. Chem. Rev. 2006, 106, 1323−1354. (12) Legrand, M. R.; Delmas, R. J. J. Geophys. Res.: Atmos. 1988, 93, 7153−7168. 12466



Article

(51) Hättig, C.; Heß, B. A. J. Chem. Phys. 1998, 108, 3863−3870. (52) Mahan, G. D. Phys. Rev. A 1980, 22, 1780−1785. (53) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420−1426. (54) Cummings, O. T.; Wick, C. D. Chem. Phys. Lett. 2013, 556, 65− 69. (55) Harrison, R. J. et al. In NWChem, A computational Chemistry Package for Parallel Computers, Version 4.1; Pacific Northwest National Laboratory: Richland, WA, 2002. (56) Kendall, R. A.; et al. Comput. Phys. Commun. 2000, 128, 260− 283. (57) Cao, J. S.; Voth, G. A. J. Chem. Phys. 1994, 101, 6168−6183. (58) Sagnella, D. E.; Tuckerman, M. E. J. Chem. Phys. 1998, 108, 2073−2083. (59) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690. (60) Dang, L. X.; Chang, T. M. J. Chem. Phys. 1997, 106, 8149−8159. (61) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577−8593. (62) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187. (63) Wick, C. D.; Dang, L. X. J. Chem. Phys. 2010, 132, 044702. (64) Wick, C. D. J. Phys. Chem. C 2009, 113, 2497−2502. (65) Wick, C. D.; Xantheas, S. S. J. Phys. Chem. B 2009, 113, 4141− 4146. (66) Roeselova, M.; Vieceli, J.; Dang, L. X.; Garrett, B. C.; Tobias, D. J. J. Am. Chem. Soc. 2004, 126, 16308−16309. (67) Roeselova, M.; Jungwirth, P.; Tobias, D. J.; Gerber, R. B. J. Phys. Chem. B 2003, 107, 12690−12699. (68) Beattie, J. K.; Djerdjev, A. M. Angew. Chem., Int. Ed. 2004, 43, 3568. (69) Wilson, M. A.; Pohorille, A.; Pratt, L. R. J. Chem. Phys. 1988, 88, 3281−3285. (70) Dang, L. X.; Chang, T. M. J. Phys. Chem. B 2002, 106, 235−238. (71) Wick, C. D.; Dang, L. X.; Jungwirth, P. J. Chem. Phys. 2006, 125, 024706. (72) Wick, C. D.; Lee, A. J.; Rick, S. W. J. Chem. Phys. 2012, 137. (73) Wick, C. D.; Cummings, O. T. Chem. Phys. Lett. 2011, 513, 161−166. (74) Wagner, C. Phys. Z. 1924, 25, 474−477. (75) Onsager, L.; Samaras, N. N. T. J. Chem. Phys. 1934, 2, 528−536. (76) Beck, T. L. Chem. Phys. Lett. 2013, 561−562, 1−13. (77) Kathmann, S. M.; Kuo, I. F. W.; Mundy, C. J.; Schenter, G. K. J. Phys. Chem. B 2011, 115, 4369−4377.

12467


HCl Accommodation, Dissociation, and Propensity for the Surface of

Recommend Documents