Aqueous NaCl and CsCl Solutions Confined in Crystalline Slit-Shaped

Dec 13, 2011 - (91, 92) For comparison, in the Supporting Information, we report the mean square displacement and the calculated self-diffusion coeffi...
0 downloads 4 Views 4MB Size
Download PDF
Article pubs.acs.org/Langmuir

Aqueous NaCl and CsCl Solutions Confined in Crystalline Slit-Shaped Silica Nanopores of Varying Degree of Protonation Tuan A. Ho,†,§ D. Argyris,†,§ D. R. Cole,‡ and A. Striolo*,† †

School of Chemical, Biological and Materials Engineering, University of Oklahoma, Norman, Oklahoma 73019, United States School of Earth Sciences, Ohio State University, Columbus, Ohio 43210, United States

‡

S Supporting Information *

ABSTRACT: All-atom molecular dynamics simulations were conducted to study the dynamics of aqueous electrolyte solutions confined in slit-shaped silica nanopores of various degrees of protonation. Five degrees of protonation were prepared by randomly removing surface hydrogen atoms from fully protonated crystalline silica surfaces. Aqueous electrolyte solutions containing NaCl or CsCl salt were simulated at ambient conditions. In all cases, the ionic concentration was 1 M. The results were quantified in terms of atomic density distributions within the pores, and the self-diffusion coefficient along the direction parallel to the pore surface. We found evidence for ion-specific properties that depend on ion−surface, water−ion, and only in some cases ion−ion correlations. The degree of protonation strongly affects the structure, distribution, and the dynamic behavior of confined water and electrolytes. Cl− ions adsorb on the surface at large degrees of protonation, and their behavior does not depend significantly on the cation type (either Na+ or Cs+ ions are present in the systems considered). The cations show significant ion-specific behavior. Na+ ions occupy different positions within the pore as the degree of protonation changes, while Cs+ ions mainly remain near the pore center at all conditions considered. For a given degree of protonation, the planar self-diffusion coefficient of Cs+ is always greater than that of Na+ ions. The results are useful for better understanding transport under confinement, including brine behavior in the subsurface, with important applications such as environmental remediation.

1. INTRODUCTION Recent advances in nanofabrication1−3 led to the development of many useful applications at the nanoscale such as nanofluidics and “lab-on-chip” processes,4−7 better understanding of biological membranes and ion-channels,8,9 and the design of ion-exclusion desalination membranes.10−12 Understanding solvent−electrolyte behavior under confinement is required to design nanoporous materials for other advanced applications, including sensors and supercapacitors. A detailed understanding of fluid−solid interactions will also provide needed insights for deploying and preserving subsurface energy systems.13 Although long simulation times are required to achieve properly equilibrated states, the growing interest in ion-exclusion processes, ion selectivity, and ion transport through pores and membranes justifies the utilization of atomistic simulations for realistic systems.14−19 Computer simulation studies have been used extensively to describe the structural properties of water near solid surfaces,20−28 and also its dynamic behavior, although to a lesser extent.29−32 The results of these investigations suggest that interfacial water properties differ significantly from those observed in the bulk. A number of experiments support these observations. In a recent review, an extensive summary is provided for both simulation and experimental investigations on this matter.33 © 2011 American Chemical Society

Regarding fundamental studies on the structure and dynamics of aqueous electrolyte solutions at interfaces, Feng et al.,34−36 Yang et al.,37 Wander and Shuford,38,39 and Shao et al.40 investigated aqueous electrolytes in carbon slit pores for developing supercapacitors. Along the same research direction, Chialvo and Cummings41 and Kalluri et al.18 investigated the partition of aqueous NaCl between bulk and slit-shaped carbon pores. Kerisit et al.42 studied the adsorption properties of monovalent ions on neutral (100) goethite surface. They pointed out that the structure of hydrated ions, the strength of ion−water interactions, and the effect of ions on water arrangement at the interface are key factors determining location and extent of ions adsorption. Argyris et al.43 simulated NaCl and CsCl solutions in fully protonated silica-based slit pores. Shirono et al.44 employed Monte Carlo and molecular dynamics simulations to study KCl distribution and transport in silica nanopores that contain both hydrophobic and hydrophilic surface patches. Their results suggest adsorption of Cl− at the silica walls and diffusion of K+ through the center of the hydrophilic pore region. Marry et al.45 suggested that the nature of counterions does not significantly alter the structure and Received: September 14, 2011 Revised: November 23, 2011 Published: December 13, 2011 1256

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

Figure 1. Side view of a partially protonated silica nanopore. Red, yellow, and white spheres represent oxygen, silicon, and hydrogen atoms of the solid substrate, respectively. The electrolyte solution contains water (red-white wireframe representation), Cl− (green spheres), and Na+ ions (purple spheres) placed between two solid surfaces.

dynamics of water near a clay surface. Bourg et al.46 studied the structure of the electrical double layer on clay surface contacting with mixed NaCl−CaCl2 electrolyte solutions. They concluded that the positions of β- and d-planes (second and third planes, respectively, formed by the adsorbate species) are independent of ionic strength or ion type. Regarding bulk electrolytes, Lee et al.47 used molecular dynamics simulations to study aqueous alkali and halides ions at ambient temperature. Diffusion coefficients48 for Na+ and Cl− ions and the solubility of NaCl and KF in water have been reported at different temperatures by Sanz and Vega.49 The effect of salt concentration on structural and transport properties of aqueous CsCl has been described by Du et al.50 These results serve as a comparison for understanding how confinement affects the ions properties. For example, the selfdiffusion coefficients for Na+, Cs+, and Cl− ions at infinite dilution in bulk water are 1.28 × 10−5, 1.88 × 10−5, and 1.77 × 10−5 (cm2 s−1), respectively, at 300 K. In our previous publication, we commented on the structure and dynamics of aqueous NaCl and CsCl solutions within fully protonated slit-shaped silica nanopores.43 The results support the existence of ion-specific effects under confinement. Because the surface was fully protonated, Cl− ions strongly adsorb onto it, while Na+ ions interpenetrate a water layer near the surface. Because of steric effects, Cs+ ions were found to accumulate in the center of the nanopore. Because of this segregation, the inplane self-diffusion coefficient for Cs+ was found to be larger than that of Na+. In the present study, we also report structural and dynamic behavior of aqueous NaCl and CsCl solutions confined within slit-shaped crystalline silica pores. Complementing our prior report, and to secure relevance with experimental observations, we interrogate the effect of the surface degree of protonation, which depends on the solution pH. The degree of protonation is considered fixed during one simulation. New advanced simulation techniques, for example, the grand canonical titration Monte Carlo ensemble 51 or the constant-pH molecular dynamics,52,53 are being developed to allow changes in protonation for titratable groups during the course of a simulation. Recent Monte Carlo simulations suggest that for pore widths comparable to those in our study, maintaining a fixed protonation state for the titratable groups is acceptable, while for narrower pores, especially in the presence of divalent cations, such an approximation may not be appropriate.54,55 For silica, the point of zero charge is pH = 2.56,57 In our simulations, the completely protonated surface is obtained by saturating all nonbridging oxygen atoms with hydrogen, yielding surface

hydroxyl groups. To mimic pH effects, we randomly remove some of the hydrogen atoms, obtaining surfaces with different hydroxyl group densities. As described in detail later, the density of hydroxyl groups on the silica substrate is controlled by a parameter, D, that varies between 0 and 1. Because simulations conducted on thin water films show that the water−structure perturbation at the interface persists for about 10−14 Å from the substrate,58 the pore width considered here (H = 12.08 Å) ensures that the confined solution is perturbed, as compared to bulk, throughout the entire pore volume. After allowing for extensive equilibration times, we calculated density profiles and self-diffusion coefficients for the ions. The results clearly demonstrate that the degree of surface protonation affects structure and dynamics of confined ions and suggest that chromatographic processes could separate NaCl from CsCl, which could have important implications for nuclear waste management and for the modification of the compositions of subsurface waters.

2. SIMULATION DETAILS Slit-shaped nanopore configurations were used in our simulations. Two silica substrates with identical surfaces were placed at a distance of H = 12.08 Å (see Figure 1) (H is the distance between the two closest layers of silicon atoms across the pore volume) along the z-axis. The (111) crystallographic face of β-cristobalite59 was used to model the solid substrate. The surface area of each periodic system is 104.8 × 100.8 Å2 (x−y plane) with a plate thickness L = 16.98 Å (see Figure 1). Details on the surface preparation can be found elsewhere.60 The size of the simulation box along the z direction, which includes the pore and the two solid substrates, is ∼4.61 nm. Periodic boundary conditions were applied along the three directions. Experimentally, the surface degree of protonation changes with solution pH.56,61 In our simulations, to study the effect of the degree of protonation on aqueous electrolyte solution properties, we adjust the hydroxyl group density according to a factor D, which ranges from 0 to 1. When D = 1, all of the nonbridging oxygen atoms are protonated, resulting in a surface hydroxyl group density of ∼4.5 OH/nm2, which corresponds to experimental data reported for amorphous silica.62 Experimentally, the surface is electrically neutral at the point of zero charge, correspondent to solution pH = 2.56,63,64 At higher pH, due to the ionization of Si−OH to Si−O−, the surface becomes negatively charged.63,64 When D = 0, all of hydrogen atoms are removed from the hydroxyl groups, and all nonbridging oxygen atoms remain exposed to the electrolyte solution. For solution 1257

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

pH above 10.5, the dissociation of the surface silanol group approaches 100%.65 Despite these changes in surface charge densities, the simulation protocol requires that the simulated system is overall neutral. To ensure that the system maintains overall electron neutrality as D changes, we adjust the charge of atoms placed 16.98 Å away from the pore surface, outside of the pore, far enough that these nonuniformly charged atoms do not interact directly with confined water. This approximation is considered analogous to immersing the system in a charged gas, as is for example done in constant-pH molecular dynamics simulations. Note that sample simulations were performed changing the thickness of the solid substrate to ensure that the results reported here are independent of the thickness of the solid substrate. For example, comparable results were obtained when the substrate thickness was 16.98 Å (present work) or 10.3 Å (previous report43). Alternatively, and maybe more realistically, we could have adjusted the number of ions in the system, but doing so we would have affected the salt concentration. We preferred not to attempt this latter strategy to compare the simulation results among systems at equal salt concentration. Sample simulations were conducted for D = 0.73 to compare the different protocols. Such comparison is provided as Supporting Information. The main result is that when the system is maintained neutral by adding Na+ ions, because of the much larger number of Na+ ions, the density peaks for Na+ are much more intense than those reported within the main text below. In particular, more Na+ ions are found near the pore center, leading to significantly higher ionic mobility than that predicted when fewer Na+ ions are present. As a consequence of the higher Na+ mobility, probably because of ion−ion correlations, Cl− ions are also found to diffuse much faster than in the simulations discussed below. These observations strengthen the main conclusion of this work; that is, the ion distribution within a narrow pore, as well as ion−water, ion−pore, and ion−ion correlations, determine their dynamic properties. Two different electrolyte solutions (NaCl or CsCl) with the same ionic strength (1 M) were confined within the silica pores. To determine the proper number of water molecules in the pore, we simulated a thin water film supported on the corresponding surface at ambient conditions and determined the number of water molecules contained within a volume equivalent to one-half of the pore volume considered in the present simulations (for details, see Supporting Information, Figure SI.1). Twice as many water molecules were then inserted in the pores. The number of water molecules changes significantly with D, reflecting changes in interfacial water structure observed for thin water films.66 As discussed in the Supporting Information, small changes in the water atomic density near the surfaces lead to significant changes in the number of water molecules that need to be simulated in the pores considered herein. The ion concentration is maintained at 1 M to ease comparison. Because of the changes in the number of water molecules, the requirement of constant ionic concentration leads to large variations in the number of ions simulated as a function of D. The number of water molecules and ions in the different systems is presented in Table 1. The CLAYFF force field67 was implemented to model the silica surface. The silicon and oxygen atoms were held at fixed positions, while the surface hydrogen atoms (for all D’s except D = 0) were allowed to vibrate and thus account for momentum exchange between aqueous solution and surface. We found this to be important when the dynamic properties of

Table 1. Number of Water Molecules and Ions in the Different Systems Simulateda surfaces D D D D D

= = = = =

0.00 0.20 0.47 0.73 1.00

pH

water molecules

ions (NaCl or CsCl)

above 10.565 3.5−10.563,64

4056 3888 3648 3190 2838

73 70 66 58 51

3.556 ∼256,63,64

a The solution pH is assumed to vary with D. The approximate solution pH for each system is also reported.

interfacial water are of interest.58 We did not account for surface reconstruction or silanol deprotonation. It is possible to directly simulate the reconstruction of silica surfaces exposed to water by implementing ab initio methods or specialized reactive force fields (not applicable to the problem of interest here because of the prohibitive computational cost associated with such methods).68−70 It should also be pointed out that realistic representations of amorphous silica surfaces could be obtained using Monte Carlo annealing.71 Such models were not implemented here because we are building on our prior work, conducted on model crystalline silica surfaces, for freestanding thin water layers and confined electrolytes. The rigid SPC/E water model was used to describe water.72 Bond lengths and angles for water were kept fixed using the SETTLE algorithm.73 This nonpolarizable water model is known to adequately reproduce experimentally observed structural and dynamic properties, for example, pair correlations and diffusion coefficients at ambient conditions.72 It has been shown that simulations of water at interfaces often require the incorporation of polarizable effects to yield realistic results.74−76 However, Sala et al.77 showed that employing either the polarizable RPOL model78 or the nonpolarizable SPC/E model for water yields results in little differences for aqueous NaCl solutions. This result, and the fact that in our system the aqueous film is confined within two silica surfaces rather than been exposed to vacuum, suggest that ignoring polarization effects should not yield pronounced errors in our results. Nonbonded interactions were modeled by means of dispersive and electrostatic forces. The van der Waals interactions were treated according to the 12-6 Lennard-Jones (LJ) potential. The LJ parameters for unlike interactions were obtained using the Lorentz−Berthelot mixing rules from pure component values.79 The electrostatic forces were described by a Coulombic potential with a cutoff set to 9 Å. Long-range interactions were calculated by the particle mesh Ewald (PME) method.80 Two-dimensional algorithms are available for treating long-ranged electrostatic interactions.81 We did not implement such algorithms because the thickness of the simulation box along the z direction, comparable to that used in our prior work, is expected to be sufficiently large to prevent nonphysical correlations among periodic replicas of the simulated system. All simulations were performed in the NVT ensemble.79 The system temperature was maintained at 300 K by using the Nosé−Hoover thermostat82,83 with a relaxation time of 100 fs. The potential parameters for sodium (Na+), cesium (Cs+), and chloride (Cl−) ions were fitted for the SPC/E water by Dang and collaborators to reproduce bulk properties.84,85 The parameters have been used by others to study ion mobility in water47 and ion adsorption at solid−liquid interfaces.86 1258

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

Figure 2. Atomic density profiles of water oxygen (a) and hydrogen atoms (b) confined within the fully protonated (D = 1), partially protonated (D = 0.73, 0.47, 0.2), and nonprotonated (D = 0) silica nanopores. Different lines represent results obtained in silica pores of different degrees of protonation D: solid black (D = 1), blue dashed (D = 0.73), red dotted (D = 0.47), green dash−dotted (D = 0.2), and purple dash−double-dotted lines (D = 0).

Figure 3. Representative simulation snapshots representing the water distribution on D = 0 (left panels), D = 0.73 (middle), and D = 1 surfaces (right). Side and top views are provided in top and bottom panels, respectively. Red, yellow, and white spheres represent surface oxygen, silicon, and hydrogen atoms, respectively. Water in different layers is shown in gray (first hydration layer), cyan, and blue (hydration layer found at ∼4 Å).

3. RESULTS AND DISCUSSION a. Water Distribution. In Figure 2a and b, we show the atomic density profiles for oxygen and hydrogen atoms of water, respectively, as a function of the distance z perpendicular to the surface. Results are obtained at different degrees of protonation D. The reference z = 0 is defined at the silicon layer of the bottom solid surface closest to the pore center. The results in Figure 2 show that the degree of protonation D dramatically affects the structure of water confined in silica pores, although the density near the pore center is comparable to bulk density found at ambient conditions. Within the D = 0 pore (purple dash−double-dotted lines), all of the surface hydrogen atoms were removed, allowing some water molecules to form the first hydration layer at z = 0.95 Å. The second and third layers were found at z = 2.15 Å and z = 4.05 Å, respectively. As D increases, more of the nonbridging O atoms in the surface are protonated. As a consequence, the position of

The equations of motion were solved using the molecular dynamics package GROMACS,87−89 by implementing the leapfrog algorithm90 with a time step of 1.0 fs. The total simulation time for all cases was 250 ns (250 × 106 time steps). Only the last 40 ns were used to calculate the properties of interest, except where indicated. We found that long simulations are necessary to obtain reliable results. Criteria for equilibration included constant energy, and constant average distribution of water and ions within the pore. The average distribution of ions in the direction perpendicular to the surface is the quantity that reaches a constant value most slowly out of the three just mentioned. When this distribution did not change over a 40 ns interval, we assumed that the system had reached equilibration. 1259

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

Figure 4. Density distributions of Na+ (purple line) and Cl− (green) ions for aqueous NaCl solutions (top panels), and Cs+ (blue) and Cl− (green) ions for aqueous CsCl solutions (bottom panels) within slit-shaped silica pores of different degrees of protonation D. The density distribution of oxygen atoms of water is also shown as gray dashed lines. Left and right y-axis represents the density of ions and water oxgyen atoms, respectively.

the first hydration layer shifts to larger z. For example, on the D = 0.47 surface, the first peak in the oxygen density profile is found at z = 1.25 Å. The positions of the second and third hydration layers do not change as D increases. Nevertheless, the intensity of these peaks decreases when D increases. As long as D is smaller than 0.73, we observe three hydration layers from the oxygen density profile. When D = 0.73 (blue dashed curve), only two well-defined hydration layers are observed, located at z = 1.65 Å and z = 4.05 Å, respectively. The position of these peaks changes as D increases further, and on the D = 1 surface the two peaks are observed at z = 2.25 Å and z = 3.9 Å. Independently of the degree of protonation D, the position of the hydration layer found at z = ∼4 Å does not change significantly, suggesting that the water structure is affected up to 4 Å from the surface. In Figure 2b, we report the atomic density profiles for water hydrogen atoms. As detailed in the case of thin water films on solid substrates,33,66 comparing the density profiles for oxygen and hydrogen atoms allows us to qualitatively assess the orientation of interfacial water molecules. For example, the atomic density profiles of water oxygen and hydrogen atoms near the surface of the D = 1 pore are similar to those observed when simulating a free-standing thin water film on a fully protonated silica surface,58 suggesting that, limited to the fully protonated case, confinement does not affect the structure of interfacial water more than a flat surface does. It is possible that as the pore width decreases this observation will no longer be valid. Data are not available for the structure of water on freestanding surfaces as a function of the degree of protonation when force fields comparable to the one used here are implemented. More quantitative assessment of the orientation of interfacial water molecules can be obtained by the calculation of appropriate order parameters, as we did for water molecules within thin films supported on various silica-based substrates. For brevity, these details are not discussed, and we instead focus on the properties of confined ions. In Figure 3, we provide simulation snapshots representing the water distribution on D = 0 (left panel), D = 0.73 (middle panel), and D = 1 surfaces (right panel). Side (top panels) and

top (bottom panels) views are provided. Red, yellow, and white spheres represent surface oxygen, silicon, and hydrogen atoms, respectively. When D = 0, three water layers are observed (see Figure 2). The water in the first hydration layer (gray) penetrates between the nonbridging oxygen atoms found on the surface, yielding a hydration layer centered at z = 0.95 Å. The second (cyan) and third (blue) layers were found at z = 2.15 Å and z = 4.05 Å, respectively. The hydration layer found at z = 0.95 Å on the D = 0 surface disappears on the D = 1 surface because the substrate is now fully protonated. The snapshots in Figure 3, right, show that the water molecules in the first hydration layer at z = 2.25 Å on the D = 1 surface (cyan) form strong hydrogen bonds (black dashed line) with the surface hydrogen atoms. On the D = 0.73 surface, the first hydration layer is found at z = 1.65 Å (cyan). Some of the water molecules in this layer seem to have features similar to those found for water molecules on the fully protonated surface (D = 1). These interfacial water molecules form hydrogen bonds (black dashed line) with the solid surface. Some other water molecules penetrate between the nonbridging oxygen atoms showing orientation similar to those of water molecules found on the completely nonprotonated surface (D = 0). On all surfaces, the snapshots confirm that the hydration layer observed at z ≈ 4 Å (blue water molecules) exhibits no change in either position or molecular orientation when D varies. b. Ion Distribution. In Figure 4 we report the density distribution of Na+ (purple line) and Cl− (green line) ions for confined NaCl solutions (top panels), and Cs+ (blue line) and Cl− (green line) ions for confined CsCl solutions (bottom panels) obtained at different degrees of protonation D. We also report the density profiles of oxygen atoms of water (gray dashed lines) for comparison. Because of the different atomic concentrations, the ionic densities are much lower than that of water oxygen. The left y-axis represents ionic density, while the right y-axis represents oxygen density. The degree of protonation affects the distribution of water molecules (described above) as well as that of the various ionic species within the pores. Note that while the density distribution of water molecules across the pore volume is always symmetric, 1260

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

when D = 0. This intense peak is less intense and wider when D = 0.2 (Figure 4b). When D = 0.47, there is evidence (small peak near the surface), suggesting that some Cl− ions move closer to the surface, and when D = 0.73 (Figure 4d) the Cl− ions clearly approach the surface. To assess by what extent the results just discussed are due to confinement, we conducted additional simulations in which a free-standing thin film of aqueous NaCl or CsCl solution was supported by a silica substrate characterized by D = 0.47. The results are reported as Supporting Information. While the positions of the various density peaks as a function of the distance from the surface were found to be consistent with the results under confinement, the peak intensities changed. c. Ionic Mobility. In Figure 5, we report the planar (x−y) self-diffusion coefficient calculated as a function of the degree of

that of the various ionic species is not always so. This is due to the nonuniform distribution of hydroxyl groups on the two pore surfaces, which in some cases lead to cooperative effects. For example, when possible, three surface hydroxyl groups act synergistically to attract Cl− ions near the surface, but in some cases there are not sufficient hydroxyl groups near each other to provide such a cooperative effect, leading to asymmetric density distributions. Considering the cations (Na+ or Cs+), the lower the degree of protonation D is, the more intense and well-defined are the peaks. For all cases, most sodium ions within the pores concentrate at z = 3.7 Å (the postion of the first oxygen peak away from the pore center), but within the D = 0 pore another peak appears at the center of the pore. The Na+ peaks within the D = 0 pore are very intense and narrow. All of the Na+ peaks become less intense, and the peak at the pore center eventually disappears as D increases. The distribution of Na+ ions within the D = 1 pore differs to some extent from what we reported earlier43 because of the slight difference in pore width (12.08 Å in this report vs 13.71 Å in our earlier work, when the pore width is calculated as the distance between the Si atoms in the surfaces across the pore volume). The Na+ ions peaks located at the center of the wider pore collapse into one peak as the pore width decreases to the pore width considered here. For the other ions considered (Cs+ and Cl−), the density distributions at D = 1 do not change significantly due to the small differences in the pore size when comparing results in ref 43 to those reported here. The distributions of Cs+ ions shown in Figure 4 are totally different when compared to those of Na+. While most of the Na+ ions accumulate not too far from the surface, all Cs+ ions remain at the center of the pore in all cases. The main change observed in the Cs+ distribution as a function of D is the intensity of the various density peaks. The peaks are found to be very intense and narrow when D = 0 and become less intense and broader when D = 1. The difference in distribution of Cs+ and Na+ ions in fully protonated pores of width 13.71 Å was discussed elsewhere.43 Despite the difference in pore width, the results in Figure 4 are in agreement with those in ref 43. They suggest that the large size of Cs+ (3.88 Å) renders it difficult to interpenetrate a water layer too close to the solid surface. Conversely, the Na+ size (2.58 Å) is comparable to that of one water molecule, and Na+ can easily be incorporated within an interfacial water layer. This observation confirms that water−ion correlations play an important role in determining ionic density distributions within narrow pores. Although the density disribution of Cl− ions depends on the degree of protonation D, our results (Figure 4) suggest that for a given D, the distribution of Cl− does not depend on the positively charged ion (compare NaCl vs CsCl solutions, top vs bottom panels, respectively). When D = 0 (Figure 4a and f), all of the negatively charged nonbridging oxygen atoms on the surface are exposed, and the Cl− ions accumulate at the center of pore because of electrostatic repulsion. When D increases, the positively charged hydrogen atoms on the surface hydroxyl groups attract Cl− ions to the surface. For example, when D = 1 (Figure 4e or j), most of the Cl− ions adsorb on the surface. The Cl− ions distribution with respect to the density distribution of water oxygen atoms in the case of fully protonated surfaces was discussed previously.43 The results in Figure 4 show a clear displacement of Cl− ions from the surface toward the pore center as D decreases from 1 to 0. As a result, an intense and narrow density peak is found at the pore center

Figure 5. Planar self-diffusion coefficients of ions and water in silica nanopores as a function of the degree of protonation D. Lines are guides to the eye.

protonation D for the ionic species as well as that for water. These data are calculated using the last 1 ns of the simulated trajectories. Note that the self-diffusion coefficients obtained under confinement are approximately 1 order of magnitude less than those typically reported in bulk liquid systems at ambient conditions. For example, the self-diffusion coefficient of bulk SPC/E liquid water at ambient conditions is (2.7−2.8) × 10−5 cm2 s−1.91,92 For comparison, in the Supporting Information, we report the mean square displacement and the calculated selfdiffusion coefficient for ions in bulk solutions. The bulk diffusion coefficients were found to be ∼1.58 × 10−5 and 1.77 × 10−5 cm2 s−1 (Na+ and Cl− in 1 M NaCl solutions, respectively), and 2.37 × 10−5 and 2.13 × 10−5 cm2 s−1 (Cs+ and Cl− in 1 M CsCl solutions, respectively). These results are in good agreement with those available in the literature, which were however obtained at infinite dilution for the electrolytes.47−50 Regarding the confined systems, the results in Figure 5 show that the degree of protonation strongly affects the in-plane selfdiffusion coefficient of confined water. The lowest self-diffusion coefficient is obtained for D = 0, and it increases as D increases. The small diffusion coefficient found in the D = 0 pore is probably related to the layered structure observed for water (see Figure 4a). As D increases, the structure of confined water becomes more disordered, which results in higher self-diffusion 1261

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

Figure 6. Planar density distributions of Na+ ions found in several layers within the slit-shaped silica pores. Panels (a) and (b) are for layers found at z = 3.75 and 5.85 Å within the D = 0 pore, respectively. Panel (c) is for the layer found at z = 3.65 Å within the D = 0.73 pore. Panels (d) and (e) are for layers found at z = 3.75 and 5.85 Å within the D = 1 pore, respectively. Density profiles in the direction perpendicular to the pore surface are shown in Figure 4 and can be used to identify the position of the layers considered here within the pores.

panel c, confirms that the ions found near the surface have reduced mobility. When D = 1, the Na+ ions show increased self-diffusion coefficients, as compared to values obtained within the other pores (see Figure 5). Within this pore type, the ions accumulate both at the pore center as well as near the surfaces. We report the planar density profiles for Na+ ions at z = 3.75 Å and z = 5.85 Å within the D = 1 pore in panels d and e of Figure 6, respectively. As expected, when the Na+ ions are close to the surface (panel d), they are very localized, leading to limited mobility, but when they are in the pore center (panel e), the density profile suggests significant mobility, which leads to enhanced self-diffusion coefficients, in agreement with our expectations. The self-diffusion coefficients obtained for Cs+ ions for a given degree of protonation D are always larger than those obtained for Na+ ions (see Figure 5). The planar self-diffusion coefficient for Cs+ is low at D = 0, consistent with the low mobility of water molecules within this pore, and it increases as D increases (blue dash line in Figure 5). For all D values considered, Cs+ ions are found to accumulate in proximity of the pore center. However, the peak intensity and width substantially change from very intense and narrow when D = 0 to less intense and wide when D = 1. The ionic distribution and mobility within a pore are due to a combination of pore−ion, water−ion, and ion−ion correlations. Our results suggest that a synergistic combination of these three tyes of correlations leads to pronounced localization and reduced mobility of the Cs+ ions within the D = 0 pore. Representative simulations were conducted for a freestanding thin film of aqueous electrolyte solution supported by the D = 0.47 surface. The results are discussed as Supporting Information, and in general show that the self-diffusion coefficients obtained under confinement are lower than those obtained in the supported film. As a control simulation, we also simulated the NaCl solution within the D = 0.73 pore when the system electroneutrality was maintained by increasing the number of Na+ ions in the system. Although qualitatively the results discussed in Figure 5 still apply for the latter system, quantitatively the self-diffusion coefficients for both Na+ and Cl− were found larger than those reported in Figure 5 because more Na+ ions were near the pore center, and because the higher Na+ mobility induced Cl− mobility via ion−ion correlations. In an attempt to explain the results discussed in Figures 5 and 6, we conducted additional characterizations of the

coefficients in the direction parallel to the surface. Because of the low ionic concentration, we expect that the water selfdiffusion coefficient does not depend on the presence of the electrolytes for the systems considered here; however, we did not consider systems without electrolytes in the present work. In our previous report,43 we concluded that the dynamic behavior of ions strongly depends on their position within the pore. For example, because Na+ accumulates closer to the surface than does Cs+ (Figure 4), the diffusion coefficient of Na+ should always be smaller than that of Cs+ within pores of the same degree of protonation. The results presented in Figure 5 are consistent with this hypothesis. The planar diffusion coefficients of Na+ ions (purple solid line) are very low at D = 0 and remain so as D increases up to D = 0.73. This suggests that the mobility of Na+ ions is reduced in all of the pores, despite changes in the mobility of water. When D = 1, the self-diffusion coefficient of Na+ ions becomes larger than that of Cl− ions, because of the different ionic distribution of the ions within the pores discussed above. To visualize the ionic mobility, in Figure 6 we present planar density distributions of Na+ ions within various planes parallel to the pore surface. The thickness of the layers considered for these calculations was 0.1 nm. These density distributions were obtained over the last 1 ns of the simulations. When the density peaks are intense and narrow, low ionic mobility is inferred. When the density peaks are less intense and broad, the results imply significant ionic mobility. The width of the highlighted areas is an indirect indication of how far the ions move during 1 ns of simulation time. In all cases considered in Figure 6, such traveled distance is rather short, which explains why the simulations used in this work need to be extremely long. The few localized red spots in Figure 6, panels a and b, suggest that the Na+ ions found in the first layer, as well as those found near the pore center, have very limited mobility within the D = 0 pore, in agreement with the low self-diffusion coefficients calculated for Na+ in this pore. Note that the mobility of water within the D = 0 pore is also very limited, which probably restricts the diffusion of electrolytes. Our results also suggest that ion−ion correlations could be responsible for the low mobility of both Cl− and Na+ ions within the D = 0 pore due to the combined effects of Cl− ion accumulation near the pore center and its overall limited mobility. Results for the planar ion−ion radial distribution function, discussed below, seem to corroborate this possibility. When D = 0.73, the Na+ ions accumulate near the surface (Figure 4d), and the planar density profile shown in Figure 6, 1262

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

Cl− (top) and for Cs+−Cl− (middle) ion pairs, respectively. For comparions, in the bottom panel of Figure 7, we report threedimensional Na+−Cl− and Cs+−Cl− radial distribution functions obtained for bulk aqueous solutions of 1 M concentration at ambient conditions. Considering the results under confinement, only the ions found within a computational volume of thickness 0.4 nm centered on the pore center were considered for the results shown in Figure 7. The results for Cs+−Cl− do not show significant changes as a function of D. It appears that some ionic accumulation takes place at center-to-center distances of ∼0.4 nm, but the change in intensity of this peak is not significant. It is instructive to point out that the most intense first peak is found for D = 0.73, in which case an approximately equal amount of Cs+ and Cl− ions are found near the pore center (see density profiles in Figure 4). The next highest peak is found for D = 0, where the density distribution for Cs+ ions is similar to that of Cl− ions (see Figure 4, bottom left panel). However, the ion−ion correlation suggested by this peak in the radial distribution function does not seem to justify the low diffusion coefficient for both Cl− and Cs+ ions within the D = 0 pore. Probably more intresting are the results for the Na+−Cl− ion pair (top panel in Figure 7). In this case, when D = 0.2, the results suggest that the ions are completely uncorrelated (the first peak in the radial distribution function is missing). This is consistent with the fact that in this pore Na+ ions adsorb near the pore surface, while Cl− ions remain closer to the pore center, and thus ion pairs cannot form easily. Significant ion−ion correlations are found in all other pore types (note the intense first peak at ∼0.4 nm and the pronounced valley between first and second peaks), but these correlations are most pronounced in the D = 0 and in the D = 1 pores. The results in the D = 1 pores are surprising because in this pore most of the Cl− ions are adsorbed onto the pore surface while Na+ ions accumulate not too far from the pore center (see Figure 4, top left panel). This observation suggests that the few Cl− ions in the pore center are strongly correlated with available Na+ ions. Should this correlation not be possible, the Cl− ions would probably adsorb on the strongly attractive surface of the D = 1 pore. The strong ion−ion correlation found within the D = 0 pore also suggests that the Na+ and the Cl− ions accumulated near the pore center (see top left panel of Figure 4), form strong ion pairs, possibly explaining the unexpected low self-diffusion coefficient calculated for both ions in the D = 0 pore. The planar diffusion coefficients of Cl− in either NaCl (green dash line in Figure 5) and CsCl (green solid line) solutions are

structure of the confined aqueous solutions. In Figure 7, we report in-plane radial distribution functions calculated for Na+−

Figure 7. In-plane radial distribution functions for Na+−Cl− (top) and Cs+−Cl− (middle) ion pairs found near the pore center. Different lines are for results obtained in pores of varying degree of protonation D. In the bottom panel, the three-dimensional Na+−Cl− and Cs+−Cl− radial distribution functions obtained in bulk 1 M aqueous systems at ambient conditions are reported for comparison purposes.

Figure 8. In-plane distribution of Cl− ions (white circles) adsorbed on the D = 1 surface (a). Planar distributions of water molecules in z = 2.25 Å (b) and z = 3.9 Å (c) hydration layers on the same surface. Planar distribution of Na+ ions within the second hydration layer on the D = 1 surface (d). See Figure 4 for the peaks position. For clarity, only a part of the simulation box is shown. 1263

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

similar, although their behavior depend significantly on the degree of protonation. The planar self-diffusion coefficient of Cl− is low at D = 0, because of the low mobility of water, and perhaps because of the association with cations, but it increases as D increases up to D = 0.73. When D increases further to D = 1, the Cl− self-diffusion coefficient decreases. The Cl− selfdiffusion coefficients do not depend significantly on the cation type because the density distributions for Cl− are for the most part indistinguishable when either aqueous NaCl and CsCl solutions are considered (Figure 4). This observation suggests that the mobility of Cl− ions is dictated primarily by pore−ion interactions, although water−ion and ion−ion correlations might also play important roles. In an attempt to quantify these correlations, in Figure 8 we report several planar density distributions obtained over the last 1 ns of the simulations. In Figure 8a, we report the in-plane distribution of Cl− ions (white circles) at contact with the surface of the D = 1 pore, shown simultaneously with the density distribution of the surface hydroxyl groups (red-green areas). It is clear that the Cl− ions are attracted to the surface because of the availability of preferential adsorption sites in which one Cl− ion simultaneously interacts with three surface hydroxyl groups. From the density profiles shown in Figure 4, it is evident that the adsorbed Cl− ions are found within the first hydration layer on the D = 1 surface, and therefore must perturb, at least locally, the structure of this hydration layer. In Figure 8b, we report the planar distribution of water in the first hydration layer. Red areas identify the preferential positions occupied by water molecules within this hydration layer. Comparing Figure 8a and b, it appears that water molecules cannot occupy the positions where the Cl− ions adsorb (these positions are highlighted by white circles), but the perturbation on the structure of the first hydration layer appears to be highly localized. In Figure 8c, we report the in-plane density distribution for water molecules found on the second hydration layer on the D = 1 surface. The density distribution appears to be patterned, in agreement with results for free-standing water films,58 although our results show that in some positions (highlighted in blue) water molecules cannot be found. Comparing Figure 8a and c, it appears that these depleted regions correspond to the sites on which Cl− ions adsorb onto the surface. Similar results, not shown for brevity, are also found for CsCl solutions in the D = 1 pore. In Figure 8d, we present the in-plane density distribution of Na+ ions found within the second hydration layer on the D = 1 surface (analyzed in Figure 8c). Comparing Figure 8c and d, it is evident that the Na+ ions do not significantly affect the structure of water within the hydration layer. It is instructive to point out that the position of the Cl− ions adsorbed on the surface (Figure 8a) does not always correspond to the position of Na+ ions within the second hydration layer (Figure 8d), suggesting that ion−ion correlations between Na+ and Cl− are only in part responsible for the distribution of these ions near the surfaces of the D = 1 pore.

ambient conditions. In all cases, the ionic concentration was 1 M. The results were quantified in terms of atomic density distributions within the pores, and the self-diffusion coefficients for ions and water along the direction parallel to the pore surface. We found evidence for ion-specific properties that depend on ion−surface, water−ion, and, only in some cases, ion−ion correlations. The degree of protonation strongly affects the structure, distribution, and dynamic behavior of water and electrolytes inside the pores. Cl− ions adsorb on the surface at higher degrees of protonation, but move to the pore center as the degree of protonation decreases. The Cl− behavior does not depend significantly on whether Na+ or Cs+ ions are present. The Na+ ions occupy preferential positions near the pore center at high degrees of protonation, and move, to some extent, toward the pore surface as the degree of protonation decreases. Because of their large size, Cs+ ions remain primarily near the pore center at all conditions considered. Because of these structural differences, the planar self-diffusion coefficient of Cs+ ions is always greater than that of Na+. When all of the nonbridging oxygen atoms on the silica surface are nonprotonated, our results show pronounced structuring of confined water, and correspondingly low mobility for water molecules and all atomic species considered. The self-diffusion coefficient for cations in the direction parallel to the surface increases as the degree of protonation increases, while that for Cl− shows a nonmonotonic behavior, with minimum selfdiffusion coefficients found at the lowest and the highest degrees of protonation considered. The results are useful for understanding and controlling ionic migration within confined geometries, with practical applications to water desalination, environmental remediation, and brine−rock interactions in subsurface systems.

4. CONCLUSIONS All-atom molecular dynamics simulations were conducted to study the structure and dynamics of aqueous electrolyte solutions confined in slit-shaped silica nanopores of various degrees of protonation. Five pore types were prepared by randomly removing surface hydrogen atoms from fully protonated crystalline silica surfaces. Aqueous electrolyte solutions containing NaCl or CsCl salt were simulated at

Financial support was provided by the U.S. Department of Energy, Office of Basic Energy Sciences, under Contract Number DE-SC0001902, and by Contract Number DE-AC0500OR22725 (Division of Chemical Sciences, Geosciences and Biosciences). Generous allocations of computing time were provided by the Oklahoma Supercomputer Center for Education and Research (OSCER) and by the National Energy Resources Supercomputer Center (NERSC). A.S. is grateful to

■

ASSOCIATED CONTENT

S Supporting Information *

Detailed explanation on the method used to determine the number of water molecules simulated within the slit-shaped pores considered here as a function of D, one sample simulation used to compare the ions distribution and mobility when the system electroneutrality is secured by changing the number of ions within the system, sample simulations for 1 M aqueous electrolyte solutions supported on one surface, and the self-diffusion coefficients for 1 M NaCl and CsCl solutions at ambient conditions. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions §

These authors contributed equally to this work.

■

1264

ACKNOWLEDGMENTS

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

(36) Qiao, R.; Feng, G.; Huang, J. S.; Sumpter, B. G.; Meunier, V. ACS Nano 2010, 4, 2382−2390. (37) Yang, L.; Garde, S. J. Chem. Phys. 2007, 126, 084706−084714. (38) Shuford, K. L.; Wander, M. C. F. J. Phys. Chem. C 2010, 114, 20539−20546. (39) Shuford, K. L.; Wander, M. C. F. J. Phys. Chem. C 2011, 115, 4904−4908. (40) Shao, Q.; Huang, L. L.; Zhou, J.; Lu, L. H.; Zhang, L. Z.; Lu, X.; Jiang, S. Y.; Gubbins, K. E.; Shen, W. F. Phys. Chem. Chem. Phys. 2008, 10, 1896−1906. (41) Chialvo, A. A.; Cummings, P. T. J. Phys. Chem. A 2011, 115, 5918−5927. (42) Kerisit, S.; Ilton, E. S.; Parker, S. C. J. Phys. Chem. B 2006, 110, 20491−20501. (43) Argyris, D.; Cole, D. R.; Striolo, A. ACS Nano 2010, 4, 2035− 2042. (44) Shirono, K.; Tatsumi, N.; Daiguji, H. J. Phys. Chem. B 2009, 113, 1041−1047. (45) Marry, V.; Rotenberg, B.; Turq, P. Phys. Chem. Chem. Phys. 2008, 10, 4802−4813. (46) Bourg, I. C.; Sposito, G. J. Colloid Interface Sci. 2011, 360, 701− 715. (47) Lee, S. H.; Rasaiah, J. C. J. Phys. Chem. 1996, 100, 1420−1425. (48) Koneshan, S.; Rasaiah, J. C. J. Chem. Phys. 2000, 113, 8125− 8137. (49) Sanz, E.; Vega, C. J. Chem. Phys. 2007, 126, 014507−014520. (50) Du, H.; Rasaiah, J. C.; Miller, J. D. J. Phys. Chem. B 2007, 111, 209−217. (51) Labbez, C.; Jönsson, B. In Applied Parallel Computing. State of the Art in Scientific Computing; Kågström, B., Elmroth, E., Dongarra, J., Wasniewski, J., Eds.; Springer: Berlin, Heidelberg, 2007; Vol. 4699, pp 66−72. (52) Donnini, S.; Tegeler, F.; Groenhof, G.; Grubmüller, H. J. Chem. Theory Comput. 2011, 7, 1962−1978. (53) Wallace, J. A.; Shen, J. K. J. Chem. Theory Comput. 2011, 7, 2617−2629. (54) Labbez, C.; Jonsson, B.; Skarba, M.; Borkovec, M. Langmuir 2009, 25, 7209−7213. (55) Barr, S. A.; Panagiotopoulos, A. Z. Langmuir 2011, 27, 8761− 8766. (56) Stumm, W. Chemistry of the Solid−Water Interface; Wiley− Interscience: New York, 1992. (57) Sverjensky, D. A. Geochim. Cosmochim. Acta 2005, 69, 225−257. (58) Ho, T. A.; Argyris, D.; Papavassiliou, D. V.; Striolo, A.; Lee, L. L.; Cole, D. R. Mol. Simul. 2011, 37, 172−195. (59) Schmahl, W. W.; Swainson, I. P.; Dove, M. T.; Graeme-Barber, A. Z. Kristallogr. 1992, 201, 125−145. (60) Puibasset, J.; Pellenq, R. J. M. J. Chem. Phys. 2003, 119, 9226− 9232. (61) Brunelle, J. P. Pure Appl. Chem. 1978, 50, 1211−1229. (62) Zhuravlev, L. T. Colloids Surf., A 2000, 173, 1−38. (63) Batteas, J.; Quan, X.; Weldon, M. Tribol. Lett. 1999, 7, 121−128. (64) Iler, R. K. The Chemistry of Silica; Wiley-Interscience: New York, 1979. (65) Philipossian, A.; Mustapha, L. J. Electrochem. Soc. 2004, 151, G456−G460. (66) Argyris, D.; Cole, D. R.; Striolo, A. Langmuir 2009, 25, 8025− 8035. (67) Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255−1266. (68) Ceresoli, D.; Bernasconi, M.; Iarlori, S.; Parrinello, M.; Tosatti, E. Phys. Rev. Lett. 2000, 84, 3887−3890. (69) Feuston, B. P.; Garofalini, S. H. J. Appl. Phys. 1990, 68, 4830− 4836. (70) Mahadevan, T. S.; Garofalini, S. H. J. Phys. Chem. B 2007, 111, 8919−8927. (71) Siboulet, B.; Coasne, B.; Dufreche, J. F.; Turq, P. J. Phys. Chem. B 2011, 115, 7881−7886.

Thanos Panagiotopoulos and Pablo Debenedetti of Princeton University, where he is currently spending his sabbatical.

■

REFERENCES

(1) Mailly, D. Eur. Phys. J. - Special Top. 2009, 172, 333−342. (2) Wu, M.-Y.; Smeets, R. M. M.; Zandbergen, M.; Ziese, U.; Krapf, D.; Batson, P. E.; Dekker, N. H.; Dekker, C.; Zandbergen, H. W. Nano Lett. 2008, 9, 479−484. (3) Quake, S. R.; Scherer, A. Science 2000, 290, 1536−1540. (4) Plecis, A.; Schoch, R. B.; Renaud, P. Nano Lett. 2005, 5, 1147− 1155. (5) Noy, A.; Park, H. G.; Fornasiero, F.; Holt, J. K.; Grigoropoulos, C. P.; Bakajin, O. Nano Today 2007, 2, 22−29. (6) Schoch, R. B.; Han, J. Y.; Renaud, P. Rev. Mod. Phys. 2008, 80, 839−883. (7) Mortensen, N. A.; Xiao, S. S.; Pedersen, J. Microfluid. Nanofluid. 2008, 4, 117−127. (8) Hille, B. Ion Channels of Excitable Membranes; Sinauer Associates Inc.: Sunderland, MA, 2001. (9) Jiang, Y.; Lee, A.; Chen, J.; Cadene, M.; Chait, B. T.; MacKinnon, R. Nature 2002, 417, 515−522. (10) Fornasiero, F.; Park, H. G.; Holt, J. K.; Stadermann, M.; Grigoropoulos, C. P.; Noy, A.; Bakajin, O. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 17250−17255. (11) Corry, B. J. Phys. Chem. B 2007, 112, 1427−1434. (12) Leung, K.; Rempe, S. B.; Lorenz, C. D. Phys. Rev. Lett. 2006, 96, 095504−4. (13) Cole, D. R.; Mamontov, E.; Rother, G. In Neutron Applications in Earth, Energy and Environmental Sciences; Liang, L., Rinaldi, R., Schober, H., Eds.; Springer: New York, 2009; pp 547−570. (14) Corry, B.; Jayatilaka, D. Biophys. J. 2008, 95, 2711−2721. (15) Kalra, A.; Garde, S.; Hummer, G. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 10175−10180. (16) Peter, C.; Hummer, G. Biophys. J. 2005, 89, 2222−2234. (17) Lorenz, C. D.; Crozier, P. S.; Anderson, J. A.; Travesset, A. J. Phys. Chem. C 2008, 112, 10222−10232. (18) Kalluri, R. K.; Konatham, D.; Striolo, A. J. Phys. Chem. C 2011, 115, 13786−13795. (19) Sint, K.; Wang, B.; Král, P. J. Am. Chem. Soc. 2008, 130, 16448− 16449. (20) Striolo, A.; Chialvo, A. A.; Cummings, P. T.; Gubbins, K. E. Langmuir 2003, 19, 8583−8591. (21) Striolo, A.; Gubbins, K. E.; Gruszkiewicz, M. S.; Cole, D. R.; Simonson, J. M.; Chialvo, A. A. Langmuir 2005, 21, 9457−9467. (22) Wang, J. W.; Kalinichev, A. G.; Kirkpatrick, R. J. Geochim. Cosmochim. Acta 2006, 70, 562−582. (23) Lopes, P. E. M.; Murashov, V.; Tazi, M.; Demchuk, E.; MacKerell, A. D. J. Phys. Chem. B 2006, 110, 2782−2792. (24) Nagy, G.; Gordillo, M. C.; Guardia, E.; Marti, J. J. Phys. Chem. B 2007, 111, 12524−12530. (25) Giovambattista, N.; Debenedetti, P. G.; Rossky, P. J. J. Phys. Chem. C 2007, 111, 1323−1332. (26) Shirono, K.; Daiguji, H. J. Phys. Chem. C 2007, 111, 7938−7946. (27) Kerisit, S.; Liu, C. X.; Ilton, E. S. Geochim. Cosmochim. Acta 2008, 72, 1481−1497. (28) Wander, M. C. F.; Clark, A. E. J. Phys. Chem. C 2008, 112, 19986−19994. (29) Gordillo, M. C.; Marti, J. J. Chem. Phys. 2002, 117, 3425−3430. (30) Striolo, A. Nano Lett. 2006, 6, 633−639. (31) Argyris, D.; Tummala, N. R.; Striolo, A.; Cole, D. R. J. Phys. Chem. C 2008, 112, 13587−13599. (32) Argyris, D.; Cole, D. R.; Striolo, A. J. Phys. Chem. C 2009, 113, 19591−19600. (33) Striolo, A. Adsorpt. Sci. Technol. 2011, 29, 211−258. (34) Qiao, R.; Feng, G. A.; Huang, J. S.; Dai, S.; Sumpter, B. G.; Meunier, V. Phys. Chem. Chem. Phys. 2011, 13, 1152−1161. (35) Qiao, R.; Feng, G. A.; Huang, J. S.; Sumpter, B. G.; Meunier, V. J. Phys. Chem. C 2010, 114, 18012−18016. 1265

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266

Langmuir

Article

(72) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269−6271. (73) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952− 962. (74) Jungwirth, P.; Tobias, D. J. Chem. Rev. 2005, 106, 1259−1281. (75) Bucher, D.; Kuyucak, S. J. Phys. Chem. B 2008, 112, 10786− 10790. (76) Picalek, J.; Minofar, B.; Kolafa, J.; Jungwirth, P. Phys. Chem. Chem. Phys. 2008, 10, 5765−5775. (77) Sala, J. J. Chem. Phys. 2010, 132, 214505. (78) Dang, L. X. J. Chem. Phys. 1992, 97, 2659−2660. (79) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 2004. (80) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (81) Yeh, I.-C.; Berkowitz, M. L. J. Chem. Phys. 1999, 110, 7935− 7942. (82) Nose, S. Mol. Phys. 1984, 52, 255−268. (83) Hoover, W. G. Phys. Rev. A 1985, 31, 1695−1697. (84) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757−3766. (85) Dang, L. X. Chem. Phys. Lett. 1994, 227, 211−214. (86) Predota, M.; Zhang, Z.; Fenter, P.; Wesolowski, D. J.; Cummings, P. T. J. Phys. Chem. B 2004, 108, 12061−12072. (87) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306−317. (88) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701−1718. (89) Hess, B.; Kutzner, C.; vanderSpoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435−447. (90) Hockney, R. W.; Goel, S. P.; Eastwood, J. W. J. Comput. Phys. 1974, 14, 148−158. (91) Berendsen, H. J. C.; van der Spoel, D.; van Maaren, P. J. J. Chem. Phys. 1998, 108, 10220−10230. (92) Mark, P.; Nilsson, L. J. Phys. Chem. A 2001, 105, 9954−9960.

1266

dx.doi.org/10.1021/la2036086 | Langmuir 2012, 28, 1256−1266