Article pubs.acs.org/JPCC
Hydronium Behavior at the Air−Water Interface with a Polarizable Multistate Empirical Valence Bond Model Collin D. Wick* Louisiana Tech University, Ruston, Louisiana 71270, United States S Supporting Information *
ABSTRACT: Molecular dynamics simulations were carried out to understand the propensity of the hydronium ion for the air−water interface with a polarizable multistate empirical valence bond (MS-EVB) model. Reasonable agreement with experiment for radial distribution functions and very good agreement for hydronium diffusion were found for the model. The polarizable MS-EVB model had no free energy minimum at the air−water interface. However, when polarizability on the hydronium ion alone was removed, a free energy of around −1.5 kcal/mol was calculated at the air−water interface. This discrepancy was found to be due to the behavior of water molecules in the first solvation shell of a hydronium ion. These water molecules contained a moderate amount of hydronium character, resulting in the delocalization of the hydronium ion. For the system with polarizable hydronium ions, this delocalization was the same at the interface as in the bulk, but for the system without polarizable hydronium ions, the delocalization increased as the hydronium approached the air−water interface. This delocalization results in a stabilization of the hydronium charge and moves it more toward the bulk, increasing its propensity for the air−water interface when hydronium ion polarizability is removed.
I. INTRODUCTION When any substance in the air reaches an aqueous droplet, it must first come in contact with an air−water interface. The interfacial structure and composition will significantly influence the accommodation of any species. For atmospheric uptake, among other things, the pH of a solution is known to significantly influence the accommodations of trace gases.1 The reason that pH can play such a role in the uptake of trace gases is that many of them are involved in acid−base chemistry. These reactions will most likely depend on not just the acidity of the droplet itself but on the acidity at its surface, which may be different than bulk. Other examples of where the acidity of the interface may be of importance is in phase transfer catalysis,2 which is known to be an interfacial reaction, and also ion complexation, in which hydronium or hydroxide ions, if they have a propensity for the interface, may displace ions to be extracted. Developing knowledge of how interfacial pH may differ from the bulk would significantly improve the understanding of the described processes and how to better control them. The topic of interfacial pH is very controversial, with differing views on if the neat air−water interface is acidic or basic. Capillary electrophoresis experiments have shown that the surface of neat water is negatively charged.3−7 Furthermore, phase-sensitive sum frequency generation (SFG) spectroscopy has found that the surface of water at an organic−water interface is apparently negatively charged, but this is different than what was found at the air−water interface.8 Newer analyses of SFG spectra have also described this result differently.9 Molecular simulation results have given a general picture in which the hydronium ion shows a free energy minimum at the air−water interface,10−15 but very recent work © 2012 American Chemical Society
has found that this depends to a degree on the strength of the water interaction with the hydronium oxygen.16 Further controversy surrounds the hydroxide anion, in which simulation and experiments have been shown to give different results depending on the method and model.11,17−24 The investigation of hydronium in solution presents a greater challenge than typical systems, due to the fact that a hydronium can share protons with adjacent waters and that a proton can transfer from one water molecule to another via the Grothuss mechanism.25,26 Because of this, classical models are generally not adequate to investigate hydronium in solution. There are different ways to model protons in solution. One of the most successful and straightforward methods is to use ab initio molecular dynamics (AIMD), which treats the system via density functional theory, and as a result proton sharing and transport are handled based on its electronic structure.15,27 The downside of AIMD is its high computational expense, usually requiring fairly small systems and short simulation times for its implementation, but larger systems have been simulated recently due to vast increases in computational power.24,28 A good compromise between computational expense and the necessary physics to investigate proton transfer events is the multistate empirical valence bond theory (MS-EVB), which was pioneered by Voth and co-workers.29−31 This method is an enhancement of the empirical valence bond method of Warshel,32 which is used to model reactions with two states. The MS-EVB model allows a proton to be shared among multiple water molecules, and in general, each state requires a Received: June 20, 2011 Revised: January 11, 2012 Published: January 17, 2012 4026
dx.doi.org/10.1021/jp209167w | J. Phys. Chem. C 2012, 116, 4026−4038
The Journal of Physical Chemistry C
Article
calculation of the total system energy. There are cross terms that are parametrized to reproduce ab initio calculations. The MS-EVB model is around an order of magnitude more computationally expensive than classical models. When dealing with ions, polarizable interactions have been found to be important in bulk water but especially when investigating the interface of water.33−39 These were predominately in relation to anions, as cations, including hydronium, have relatively low polarizabilities. Moreover, the addition of polarizable interactions increases the computational expense of simulations, even though there are ways to substantially mitigate this.40,41 Polarizability is also known to affect the dynamics of molecular systems42 and has been used to produce an MS-EVB model that provides very good agreement for hydronium diffusion in water.43 While polarizability has not shown a significant influence on alkali cation interfacial behavior, these cations have always been repelled from the air−water interface.37,44 A cation such as hydronium, that has a propensity for the air−water interface, may be influenced to a degree by polarizability. An air−water interface does have a symmetry breaking between a polar and a nonpolar medium. This symmetry breaking would be expected to result in a change in the average induced dipole for an interfacial species with respect to the bulk. This may have an influence on interfacial behavior. The only example, to my knowledge, that calculated the free energy difference between the air−water interface and water bulk with a model that utilized both polarizability and proton sharing was the work by Lee and Tuckerman in an AIMD investigation.15 Their work showed an interfacial free energy slightly more negative than the bulk free energy. This present paper details the development of a polarizable MS-EVB hydronium ion model and its implementation for the determination of the influence of hydronium polarizability on its interfacial free energy.
Table 1. Potential Parameters for the Flexible Water Model and Hydronium Anion Used in the MS-EVB Potential H2O value
H3O+ value
D0 (kcal/mol) α (Å−1) r0 (Å) ka (kcal mol−1 rad−2 Å−1) kb (kcal mol−1 rad−2 Å−2) kc (kcal mol−1 rad−2) θa (rad Å−1) θb (rad Å−2) θc (rad) σO (Å) εO (kcal/mol) σH (Å) εH (kcal/mol) polarizability qH (e) qO (e)
105.69 2.279 eq 2 51.91 36.8 29.3 0.927 −1.361 1.303 3.67 0.36 0.0 0.0 1.444 eq 6 eq 6
144.12 1.85 0.98 0.0 0.0 89.15 0.0 0.0 1.9635 3.56 0.36 0.0 0.0 0.98 0.4722 −0.4166
interactions was due to the strained geometries often found in reactive configurations. k UBend = 0 (θ − θ0)2 2
(3)
k 0 = ka(r1 − r2) + kbr1r2 + kc
(4)
θ0 = θa(r1 − r2) + θbr1r2 + θc
(5)
where the constants, ka and θa, along with b and c parameters for these are given in Table 1. These were parametrized to reproduce ab initio calculations 45 and compared with experimental IR spectroscopy, giving peaks at 1595, 3610, and 3710 cm−1, which agree reasonably well with the experimental values of 1595, 3657, and 3756 cm−1.46 A further parametrization was carried out to determine how the charge distribution of the water atoms depends on its intramolecular geometry.47 The total water charge was kept neutral, but depending on the HOH angle and the OH distances the hydrogen charge was able to change. The hydrogen charge is determined by
II. MOLECULAR MODEL DETAILS II.A. Flexible Water Model. A water model that closely mimics a previously developed model was used for this work.21 The model is flexible with four sites and includes polarizability. There are three atomic sites, two hydrogen sites, one oxygen site, and an m-site located along the bisector of the oxygen− hydrogen bonds. For bond stretching, a Morse potential is used, in which the two oxygen−hydrogen bond stretches are coupled.
qH = − 0.14rOH + 0.14θ HOH + 0.398
(6)
where θ is in radians, and the opposing charge is located on the m-site (not the oxygen atomic position), which adjusts so the water molecule charge is neutral. At the equilibrium gas-phase geometry (r = 0.9572 Å and θ = 104.52°),48 the hydrogen charge is 0.519e, giving the correct gas-phase dipole of 1.85 D. In the liquid phase, the average bond length is 0.97 Å, and the average bend angle is 106.1°, in good agreement with experiment.49 The location of the m-site is dependent on the geometry of the water molecules and is given by the following equation
UBond = D0[1 − e−α(r1− r0,1)]2 + D0[1 − e−α(r2 − r0,2)]2 (1)
where r0,1 and r0,2 are given by r0,1 = 0.91656 + 0.01159 × r2 , r0,2 = 0.91656 + 0.01159 × r1
property
(2)
xm = γ(x H1 + x H2) + (1 − 2γ)xO
and r1 and r2 represent the distances between each hydrogen and the oxygen. The equilibrium bond length of r1 is dependent on the value of r2, and vice versa. Table 1 also gives the values for α and D0, which were fit to reproduce the symmetric and antisymmetric stretches for the gas-phase IR spectroscopy of water. The hydrogen−oxygen−hydrogen bond bending potential is dependent on the water bond lengths. The reason for using a somewhat complicated procedure for the bonded
(7)
where the value for γ was set to 0.182. As stated, this water model is quite similar to a previously developed one,21 but this model utilizes a different interaction for atomic repulsion and van der Waals interactions. While the previous model utilized a Lennard-Jones potential, the current one uses a Buckingham exponential-6 (exp-6) potential, in which a single site is located on the oxygen atomic position. It should be noted that the exp4027
dx.doi.org/10.1021/jp209167w | J. Phys. Chem. C 2012, 116, 4026−4038
The Journal of Physical Chemistry C
Article
Figure 1. Schematic of the implementation of the MS-EVB model for a hydronium cluster.
⎯r → ⎛ ⎞ 3→⎯ ij rij 1⎜ Tij = 3 f2 (rij) 2 − f1 (rij)⎟ , ⎟ rij ⎜⎝ rij ⎠
6 potential greatly improves agreement with experiment for the water structure (as discussed in Section IV). The exp-6 potential is given by ⎡ 6⎤ ⎛ ⎡ 6 r ⎤⎞ ⎛ σ ⎞ Uvdw = ε⎢ exp⎜λ⎢1 − ⎥⎟ − ⎜ ⎟ ⎥ ⎝ ⎣ σ ⎦⎠ ⎝ r ⎠ ⎥⎦ ⎢⎣ λ
⎡ q q f (rij) ⎤ i j0 ⎥ − 1 ∑ μ⃗ E 0 i i 4πε0 ⎥⎦ 2 i
