J. Phys. Chem. C 2009, 113, 2497–2502
2497
NaCl Dissociation Dynamics at the Air-Water Interface Collin D. Wick Department of Chemistry, Louisiana Tech UniVersity, Ruston, Louisiana 71270 ReceiVed: September 5, 2008; ReVised Manuscript ReceiVed: December 15, 2008
The rate for NaCl dissociation was compared between the air-water interface and water bulk using the transition path sampling formulism to begin understanding how an interface affects reaction dynamics. The results showed that the dissociation dynamics at the interface was much slower than in the bulk. Free energy of dissociation calculations were carried out utilizing umbrella sampling for NaCl dissociation at the interface, and a higher barrier to dissociation was observed, which is consistent with the slower interfacial dissociation dynamics. Using the dissociation free energy profile at the interface, transition-state theory was used to predict a rate constant that was lower than the one calculated using transition path sampling. This led to the conclusion that it is difficult to capture the complex nature of interfacial reactions with mean free energy profiles. Analysis of transition states harvested from the simulations found that the NaCl vector was significantly less aligned to the normal of the air-water interface at the transition state than on average and that the Cl- ion had a high probability to be found at the Gibbs dividing surface at the transition state. I. Introduction Interfaces are a barrier to molecule transfer and often serve as a distinct environment for reactions that affect atmospheric composition, biological processes including drug delivery, liquid chromatography, and waste cleanup, to name a few. Waste cleanup and ion complexation are likely dependent on interfacial complexation reactions, and interfacial reactions strongly influence atmospheric reactions in aqueous aerosols. Understanding the microscopic structure of interfaces and how they affect chemical reactions will be of great benefit to these and countless applications. Probing dynamics at the interface can be a great challenge. Molecular transport and dynamical reaction properties at interfaces have been investigated by second harmonic generation spectroscopy,1,2 and, to the knowledge of the author, an interfacial reaction rate constant has only been determined experimentally for the Cl2- radical attack of ethanol at the air-liquid interface.3 Also, molecular beam studies can bring insight into interfacial reaction kinetics.4 However, a molecular level understanding of many interfacial reaction mechanisms and dynamics, and how they affect one another, are far from apparent. Computational methods can provide a pathway for a strong understanding of processes on the molecular level, and also can be used calculate dynamical properties. The computational determination of reaction rates can be achieved by a variety of techniques, including variational transition-state theory,5 reactive flux formulism,6,7 and Grote-Hynes theory,8 to name a few. However, in the case of many interfacial reactions, the reaction mechanism is not known a priori, including the dividing surface, making traditional methods for calculating reaction rate constants difficult to implement. Having an interface present adds greater dimensionality to a system in which its effects on reaction mechanism are not known a priori. Previous computational work to understand the role of an interface on reactions have investigated a model SN1 reaction at the air-water and organic-water interfaces,9,10 in which a significant increase in barrier height was observed at the interface. Also, ICN photodissociation was investigated, finding that the air-water interface significantly affects its reaction dynamics and vibrational
relaxation due to the lower water density at the interface.11 For interfacial NaCl dissociation at an aqueous-organic interface, it has been investigated using the potential of mean force technique for nonpolarizable models.12 However, when using nonpolarizable models, ions are generally repelled from the air-water interface.13 When polarizability is included, it has been found, though, that many anions, depending on size and polarizability, can have a propensity for the air-water interface.14-20 The current study will investigate a model reaction at the air-water interface: the dissociation rate constant of NaCl. In bulk, NaCl dissociation has been studied using transition path sampling (TPS),7,21 and gave similar results as Grote-Hynes theory, and identical results as reactive flux formalism.7 In addition, the rate constant for the transport of a hydroxyl radical from the air-water interface into water bulk has been calculated using TPS, and gave similar results as the Grote-Hynes theory.22 For NaCl dissociation, the interatomic NaCl distance can be chosen as a reaction coordinate, and for hydroxyl radical transfer into the bulk, the normal to the interface can be chosen. However, for NaCl dissociation at the interface, the extra dimensionality the interface adds makes a determination of the reaction coordinate much more difficult. Specifically, the NaCl interatomic distance and the distances from the Gibbs dividing surface (GDS) for each of the ions would need to be incorporated into the reaction coordinate, as they may all influence the dissociation dynamics. This will hinder the use of Grote-Hynes theory, which requires the identification of a free energy maximum along the reaction coordinate. In contrast, the TPS technique does not require knowledge of the reaction coordinate but only of well-defined reactant and product states,23 after which it can be used to calculate a rate constant.24 Because of this, TPS is an attractive method for the calculation of reaction rate constants when a reaction coordinate is difficult to define. This article outlines the investigation of NaCl dissociation at the air-water interface using the TPS formulism to calculate a reaction rate. In addition, insight is given into the interfacial energetic and structural properties that influence the interfacial NaCl dissociation mechanism and rate.
10.1021/jp807901j CCC: $40.75 2009 American Chemical Society Published on Web 01/21/2009
2498 J. Phys. Chem. C, Vol. 113, No. 6, 2009 II. Simulation Methodology A. Molecular Models. To model water, the rigid and polarizable Dang-Chang force field was utilized, which includes a single Lennard-Jones site at the oxygen positions, two positive charges at the hydrogen atomic positions, and a negative charge and point polarizability on an m-site, located along the oxygen hydrogen bisector.25 For the simulation of NaCl, a polarizable model was used also.26 Because the described simulations were developed to compare NaCl dissociation at the interface and bulk, polarizability is of utmost importance for the models, as anion interfacial activity is significantly affected by the inclusion of polarizability.26 B. TPS Formalism. The TPS formalism has been described in detail previously,23 and only a brief account will be given here. TPS requires parameters for reactant and product states to be defined and for dynamical trajectories to be limited to only those that connect these reactant and product states. An initial trajectory is required that connects the reactant and product states (which is described in section II.C. for this work), and the TPS formulism induces the creation of additional novel trajectories connecting reactant to product. The creation of new trajectories is carried out by performing random transitions from one trajectory to the next by either the shooting or shifting moves.23 The shooting move randomly picks a position in the trajectory, randomly changes the momentum of the Na+ and Cl- ions based on a Boltzmann distribution (at 298 K for this work), and runs trajectories in the forward and reverse directions based on the new momenta. If the new trajectory connects the reactant and product states, then it is accepted; otherwise it is rejected and another move is performed. The shifting move simply shifts the trajectory either forward and backward in time by a random amount and is accepted if the shifted trajectory connects reactant and product states. One major key for the implementation of the TPS formalism is the definition of reactant and product states. For NaCl dissociation, the interatomic distance showed a free energy barrier at 3.7 Å in bulk.27 Because of this, the reactant state was chosen to be less than 3.2 Å, and the product state was chosen to be greater than 4.6 Å. For dissociation at the interface, two questions have to be resolved. The first is if the dissociation will have a free energy barrier at a similar position. The effect of the interface has been investigated for the free energy of NaCl dissociation previously, and the barrier position with respect to NaCl interatomic distance was found to be similar at the interface and in the bulk.12 Because of this, the NaCl interatomic distance for the reactant and product states was kept the same at the interface as in the bulk. The second issue is how to incorporate the position of the interface itself with an interfacial reaction. To define an interfacial reaction, the reactant state was required to have Cl- within 3 Å of the GDS of water. This was chosen because Cl- has an enhanced density with respect to bulk water in this region for a 1 M NaCl aqueous solution.28 No constraint on the reactant state was imposed as far as the Na+ position with respect to the GDS, and no constraints on the product states for the ion position with respect to the GDS was imposed. The reason for choosing only the Cl- interfacial coordinate and not the Na+ position is that constraining the Na+ to be near the interface would result in starting configurations that are not typical (because it is repelled from the air-water interface).28 In summary, the reactant state for an interfacial reaction was defined by NaCl interatomic position and the Cl- position with respect to the GDS, and the product state was only defined with respect to the NaCl interatomic position. It is important to test
Wick these assumptions to make sure the proper reactant and product states are defined, and this is described later. To determine the rate constant of a reaction, the timederivative of the autocorrelation function C(t) for the creation of products from reactants is calculated,
C(t) )
〈hA[rNaCl(0)]hB[rNaCl(t)]〉 〈hA[rNaCl(0)]〉
(1)
where hA represents the characteristic function for the reactant state, and hB is for the product state. These functions are 1 when true at the specified time and 0 otherwise. For the calculation of the rate, the time derivative of C(t) can be factorized into two terms,24
k(t) )
dC(t) 〈dhB ⁄ dt〉 ) × C(t′) dt 〈hB(t′)〉
(2)
where t′ is a fixed time chosen usually to be much less than the total time of the trajectory sampled, and hB is defined as the characteristic function equaling 1 if a trajectory originating in the reactant state is in the product state in the interval [0, t], and 0 otherwise. For the evaluation of hB, trajectories that begin in the reactant state and Visit the product state at least once during the trajectory are created. For this work, a total of 20 ,000 trajectories with lengths of 2 ps for the bulk trajectories and 4 ps for the interfacial trajectories were carried out for each set of the systems to evaluate hB(t). A second set of trajectories are required to calculate C(t′). This was achieved by calculating the probability distribution with respect to reaction coordinate (which is chosen as the NaCl interatomic distance) for a set time (t′). This time is usually chosen to be much less than for the first set (used for the calculation of hB(t)) trajectory time length, because more trajectories are usually required for its evaluation, and was chosen to be 0.4 ps for this work. The probability distribution is then integrated over all regions in the product state (for NaCl interatomic distances greater than 4.6 Å) to determine C(t′). To improve the sampling of this distribution, umbrella sampling was carried out in which the trajectories were constrained to finish in different NaCl interatomic distance windows.24,29 Using the TPS technique, the trajectories were restricted to all begin in reactant state and to finish in the defined windows. Two separate windows were chosen, one for NaCl distances less than 4.2 Å, and one for distances greater than 4.0 Å. A total of 16 000 trajectories were carried out for each window for both the interfacial and bulk systems. The probability distributions of the adjacent windows are matched when they overlap to create a probability distribution over the entire trajectory. C. Simulation Details. Two different types of systems were simulated. One system was in a cubic box, in which 440 water molecules were simulated. The other had 550 water molecules, and was in a slab geometry, in which a rectangular box was elongated in the z direction, with two vapor-liquid interfaces formed bisecting the z axis (or the z axis was normal to the interfaces). The dimensions for the elongated box were 25.5-25.5-75.0 Å for the x-y-z axes, respectively, with the liquid occupying approximately 26.6 Å along the z axis from interface to interface, and the liquid and vapor were periodic in the x and y axis. Both simulations were initially equilibrated at a constant temperature of 298 K using Nose-Hoover chains,30 and the cubic box had its pressure set to 1 atm, which was set using the Berendsen barostat.31 All dynamical calculations were carried out in the NVE ensemble. The water geometry was kept rigid using the SHAKE and RATTLE algorithms.32 A potential truncation of 9 Å was enforced for LJ interactions, and the
NaCl Dissociation Dynamics
Figure 1. Probability distribution with respect to NaCl distance at the interface and in the bulk for a 0.4 ps trajectory. The vertical lines represent the barriers for reactant, intermediate, and product states.
particle mesh Ewald summation technique was used for longranged electrostatics.33 The generation of the first trajectory to be used by the TPS technique proceeded by simply placing a NaCl pair together at the interface (or in the bulk for the bulk simulations), running molecular dynamics simulations, and waiting until a dissociation trajectory occurred due to random fluctuations in the allotted time (6 ps for the interface and 2 ps for the bulk). III. Results and Discussion A. Probability Distribution for Short Trajectories. The probability distribution for 0.4 ps trajectories to be at certain NaCl interatomic distances for the bulk and interfacial regions are given in Figure 1. For reference, the reactant state has an interatomic distance less than 3.2 Å, and the product state is at an interatomic distance of greater than 4.6 Å. Both distributions have qualitatively similar shapes, with a region of high probability around 2.7 Å, followed by a decrease, then a flattening, and finally a monotonic decrease. The difference between the two distributions is mainly in the fact that the distribution for the interfacial reaction drops off much more rapidly even before the intermediate state (or an interatomic distance of 3.2 Å) is reached. This clearly shows that it is more difficult for the intermediate state to be reached at the interface. Furthermore, because of this, the product state has a much lower probability of being reached at the interface than in the bulk, which will reduce the reaction rate for the interface. It should be noted that this analysis is for only very short trajectories (0.4 ps) and is only one factor used for the rate determination (C(t′)) in eq 2. B. Rate Constant Determination. The rate constant as a function of time is given in Figure 2 for the bulk and interface. The bulk rate constant plateaus relatively quickly at around 1.5 ps, giving a rate constant of 0.11 ( 0.01 ps-1, which compares to the rate constant for nonpolarizable NaCl dissociation of 0.05 ps-1,7 showing a faster rate for polarizable ion dissociation. It should be noted that there are multiple differences between these models, including the water model used, and that polarizability may not be the only factor that is responsible for the increased rate. Previous simulations of polarizable NaCl dissociation (but with a different water model than the current work) using Grote-Hynes theory found a rate constant of 0.087 ps-1,27 which is very similar to the results given here. Whereas there is debate as to the quality of Grote-Hynes theory in predicting rate constants for NaCl dissociation,7 it appears that it gives reasonable results for NaCl dissociation, as well as estimating the hydroxyl radical residence time at the interface.22 Grote-
J. Phys. Chem. C, Vol. 113, No. 6, 2009 2499
Figure 2. Rate with respect to time, with the plateau values represented by the horizontal lines.
TABLE 1: Average Values of Selected Properties in Reactant, Product and Intermediate States from Trajectories Used to Determine hB(t)a 〈µ〉Cl (D) 〈µz/µ〉 Cl-b 〈zNa+〉 (Å)c 〈zCl-〉 (Å)c -
reactant state
product state
intermediate state
2.21 0.774 3.21 2.01
2.21 0.784 3.62 1.91
2.11 0.666 4.12 2.61
a
The subscripts represent the uncertainty in the last digit. Represents the cosine of the angle of the induced dipole with the z axis (positive values represent a dipole pointing toward the liquid). c Average distance from the GDS of water. b
Hynes theory, though, would be much more challenging to use for an interfacial reaction because there is no well-defined transition state due to the added dimensionality of the reaction trajectory from the presence of the interface. At the interface, the rate constant takes a much longer time to plateau than in the bulk. This suggests that NaCl dissociation has a much longer memory effect at the interface than in the bulk. This is probably due to the fact that NaCl dissociation is much slower at the interface. Also, because of the added dimensionality of the interface, a greater number of different trajectories are possible at the interface because only the reactant state has an additional constraint on the anion position. For instance, Na+ may move perpendicularly with respect to the interface into the bulk, or more parallel to the interface. In addition, Cl- may transfer into the bulk, or stay at the interface itself. These factors may give rise to very diverse trajectories (which will be described in detail below), each with different rates, and to determine the overall rate of the NaCl dissociation at the interface requires sampling nearly all of them. After the plateau is reached, the rate constant for NaCl interfacial dissociation is 0.026 ( 0.005 ps-1, which is less than a quarter as high as the rate constant in the bulk. This is consistent with what was expected from the NaCl distance probability distributions for a 0.4 ps trajectory because it was found that there was a much lower probability for an interfacial NaCl trajectory to reach the product state. C. Average Reaction Properties. To understand the pathways for interfacial NaCl dissociation trajectories, Table 1 gives a variety of properties calculated for the reactant, product, and intermediate states from trajectories used to calculate hB(t). It should be noted that the intermediate state values are not transition states but simply cases where the NaCl interatomic distance is within the intermediate region (i.e., between 3.2 and 4.6 Å). The Cl- ion has a dipole that is pointing toward the liquid bulk to the same degree for the reactant and product states. In addition, Cl- has a similar z-position with respect to the GDS
2500 J. Phys. Chem. C, Vol. 113, No. 6, 2009
Wick
Figure 3. Free energy as a function of NaCl distance with the Clnear the interface with its position near 2.0 Å and 0.0 Å from the GDS surface.
of water for both states. However, in the intermediate state, the Cl- ion appears to migrate more toward the liquid interior by around 0.5 Å, and its dipole points slightly less toward the interior than in the product and reactant states. The position of the Na+ ion follows similar trends as the Cl- ion in the intermediate state, in that it migrates toward the interior when compared with the product and reactant states. This gives a picture of an interfacial NaCl pair migrating somewhat toward the interior as it dissociates. However, this does not include specific information about the transition state, which may be significantly different. D. Interfacial Free Energy Profile. To compare the free energy profile of the interface with the bulk, umbrella sampling29 was carried out to understand the free energy with respect to NaCl distance when Cl- is near the interface at 298 K. To achieve this, two harmonic umbrella potentials were utilized, one for the NaCl interatomic distance (rNaCl), and one for the Cl- z-position,
U ) k0(r - r0)2 -1
-2
(3)
where k0 was 5.0 kcal mol Å for both rNaCl and zCl, and two situations were simulated with different zCl r0 values. One case was simulated in which r0 ) -3.0 Å for the zCl value, with respect to the GDS, but only for cases in which zCl was less than -3.0 Å, so all cases where the Cl- position was closer to the GDS than -3 Å had the same umbrella energy with respect to zCl. For the second case, a value of r0 ) 0 was used for zCl. With greater than 99% probability, the Cl- position was within 1 Å of the umbrella limits or a zCl value greater than -4 Å for r0)-3.0 Å and between -1 and 1 Å for r0 ) 0. For rNaCl, r0 was varied, with six separate 800 ps simulations carried out, each with a different r0 values, including 2.5, 3.0, 3.5, 4.0, 4.5, and 5.0 Å. The probability distributions for each simulation (corrected for the fact that an umbrella sampling potential was used) were matched in overlapping regions to make probability curves that stretched from around 2.3 Å to 5.5 Å, in which a free energy profile (W(r)) was created, and is given in Figure 3. The free energy maximum is at an interatomic distance of 3.5 to 4.0 Å, in agreement with what was found in the bulk.27 Also, this shows that the boundaries for the NaCl dissociation reactant and product states at the interface should be similar at the interface and bulk. In contrast, the free energy barrier is approximately 2.8 kcal/mol (with the energy values having an uncertainty of 0.2 kcal/mol by comparing four independent simulations) for z0
