J. Phys. Chem. B 2009, 113, 13993–13997
13993
Computational Studies of Aqueous Interfaces of SrCl2 Salt Solutions Xiuquan Sun,† Collin D. Wick,‡ and Liem X. Dang*,† Chemical and Materials Sciences DiVision, Pacific Northwest National Laboratory, Richland, Washington 99352, and Chemistry Department, Louisiana Tech UniVersity, Ruston, Louisiana 71270 ReceiVed: July 1, 2009
The electron density profiles and corresponding surface structures of an aqueous interface of SrCl2 salt solution were computed by use of molecular dynamics simulations. We used both polarizable and nonpolarizable potential models to describe molecular interactions. The results demonstrate that the polarizable models captured the essential features of the corresponding X-ray reflectivity experimental data while the corresponding nonpolarizable models could not. In addition, we demonstrated that the shape of the X-ray reflectivity curve could be quantitatively reproduced if the simulations were carried out at lower SrCl2 concentrations, making it likely that the polarizable models used in this study somewhat overestimate the surface concentration of ions. However, significant interfacial enhancement of both Sr2+ and Cl- appears necessary to reproduce the experimental spectra. This is in contrast to systems with monovalent cations, which have generally been found to have a double layer, in which anions are enhanced at the surface but cations are repelled. I. Introduction Ion distributions at aqueous/vapor interfaces have attracted intensive attention due to their fundamental and practical importance. The traditional view is that the liquid/vapor interface is devoid of ions based on interpretations of surface tension measurements1–4 due to the Gibbs adsorption isotherm.5–8 This argument is challenged from both theoretical and experimental studies, with an enhancement of anion concentrations at the liquid/vapor interface being observed from both molecular dynamics (MD) simulations9–19 and spectroscopy.20–26 One recent effort on this topic is the work of Sloutskin et al.,27 who used X-ray reflectivity techniques to explore the surface structure of concentrated aqueous salt solutions. Their results are a very important contribution to the understanding of liquid/vapor interfaces, and if they can be reliably interpreted, they can bring major improvements to the understanding of interfacial ion behavior. Molecular modeling approaches, such as molecular dynamics simulations, are powerful tools with the ability to probe the ion behavior at liquid interfaces, providing a strong microscopic understanding of their structures, distributions, and behavior.28 Unfortunately, observations made from MD simulations often cannot be directly compared with experimental measurements, making conclusions extracted from them difficult to verify. Because of this, effort has been put forth to bridge the gap between experimental measurements and MD simulations for aqueous salt solutions. For instance, sum frequency generation spectra for salt solutions can be generated from MD simulations,29,30 but these are somewhat indirect probes of ion distributions, as they primarily describe water vibrations. X-ray reflectivity spectroscopy, though, probes the structure of interfacial electron density, and recent work has been carried out to compare MD simulation results with spectroscopic measurements by Sloutskin et al.27 However, the MD simulations were not performed at the same concentrations for RbBr or with the * Corresponding author. † Pacific Northwest National Laboratory. ‡ Louisiana Tech University.
same cations for SrCl2 as experiment, requiring the simulation results to be scaled or approximated to correspond with the experiments. Consequently, the outcome of such comparisons cannot be deemed rigorous. Very recently, a direct comparison with X-ray reflectivity experiments using molecular dynamics simulations for concentrated RbBr aqueous solution has been carried out.19 The ion concentrations were found to be enhanced at the interfaces, and the calculated structure factor reproduced the experimental one qualitatively from simulations with polarizable potential models, while the simulations with nonpolarizable potential models showed contrasting results. In this paper, we describe our investigation of the distribution of the SrCl2 and NaCl ions at the aqueous/vapor interfaces using MD simulations. To the best of our knowledge, this work is the first MD calculation of the X-ray reflectivity of the concentrated divalent salt solutions. In addition to the potential models, divalent cations behave significantly differently than monovalent ones due to their strong ion-ion interactions. As a consequence, many body effects may have enhanced importance for these systems. One of our main goals is to make a direct comparison to the measured data, and we also discuss ways to improve the agreement between MD and X-ray reflectivity, including new interpretations of the aqueous salt surfaces. We used both polarizable and nonpolarizable potential models and full 3D periodic Ewald summation to investigate the distribution of ions at the aqueous interfaces. II. Computational Methods Both nonpolarizable and polarizable potentials were used to describe the molecular interactions in our simulations. The nonpolarizable potentials consisted of electrostatic and van der Waals interactions, with the TIP4P31 water model used for our simulations. For the polarizable simulations, the DC9732 water model was used. Atomic polarizabilities were present on all salt ions for the polarizable simulations, and all of the ion parameters are given in previous work.33,34 The systems in our simulations consisted of 4000 H2O molecules, 180 Sr2+ and 360 Cl- ions, corresponding to a 2.7 M SrCl2 aqueous solution under ambient
10.1021/jp9079525 CCC: $40.75 2009 American Chemical Society Published on Web 09/29/2009
13994
J. Phys. Chem. B, Vol. 113, No. 42, 2009
Sun et al.
conditions. The simulations were carried out in boxes with dimensions of 52 × 52 × 150 Å3, with liquid occupying approximately 50 Å of the simulation box along the z-axis and vapor occupying the remainder of the box, resulting in two liquid/vapor interfaces. Analytical tail corrections were enforced for the van der Waals interactions, and the particle mesh Ewald summation technique was used to evaluate long-range electrostatic interactions.35 We employed the Berendsen heat bath coupling of 2.5 ps to keep the temperature constant near 300 K.36 After a 10 ns equilibration period with a time step of 2 fs, the data for analysis were averaged over 10 ns in the NVT ensemble. It should be noted that we found increases in surface tensions with respect to neat water for both nonpolarizable and polarizable simulations of 2.7 M SrCl2, in agreement with expected behavior.1–4 The detailed description of the calculation of X-ray reflectivity can be found in other papers.19,27,37 Here, we will only summarize some concepts. The definition of X-ray reflectivity is given by
R(qz) ) RF(qz)|Φ(qz)| 2exp(-σe2qz2)
(1)
where qz is the surface-normal wave vector, RF(qz) is the Fresnel reflectivity (i.e., the reflectivity of a flat and continuum surface), and Φ(qz) is the structure factor. σe is the effective amplitude induced by thermal capillary wave and nonthermal zero-point motion, including the shape of the molecule. The structure factor Φ(qz) can be calculated from the electron density profile F(z) along the surface normal direction z: ∞
Φ(qz) )
∫[
dF(z) 1 exp(iqzz) dz Fb -∞ dz
]
(2)
and Fb is the bulk electron density. The total electron density was computed from the sum over the product of number densities ni(z) along surface normal direction, and atomic number Zi of species i, F(z) ) ∑i ni(z)Zi. III. Results and Discussion Figure 1 shows the normalized number density profiles of H2O, Sr2+, and Cl- calculated from the simulations by use of polarizable potentials. Two interfaces are located at 50 and 100 Å along the z direction. The structures of the interfaces are symmetric, which gives evidence for the stability for the structural properties of the interfaces. The most striking result shown by Figure 1 is that the ions are highly enhanced at the liquid/vapor interfaces, with the number densities of both ions at the surface being roughly 5 and 8 times those of the bulk for Cl- and Sr2+, respectively. These results are somewhat unexpected since the solvation energy of the Sr2+ ion is very strong (≈-350 kcal/mol)38 and thus one would expect that the Sr2+ and Cl- ions should be localized in the bulk instead of at the interface. The peak shape and positions of the ions are almost the same except that Sr2+ has a slightly narrower peak than Cl-. From the MD simulations with polarizable potentials, the anions in monovalent salt aqueous solutions are enhanced at the liquid/vapor interface while the cations are depleted, but the ratio is reversed a few angstroms toward the liquid interior, causing a double layer to form.9–16 In this study, divalent Sr2+ cations have similar enhancement at the surface as Cl-, in contrast to the double layers formed with monovalent cations. Figure 1 also gives a snapshot of a simulation with polarizable
Figure 1. (Top) Snapshot from the simulations with the polarizable potential. Water molecules are in blue. The Cl- ions are in green and the corresponding Sr2+ ions are in red. (Bottom) Computed normalized number density profiles for the 2.7 M SrCl2 aqueous salt solution interface by use of polarizable potential models for water, Sr2+ and Cl-.
potentials, showing a very large number of ions present at the surface. Unlike the monovalent ions, which form solventseparated ion pairs in dilute solutions, Sr2+ and Cl- primarily exist as contact ion pairs in the solution due to the large net charge on Sr2+. In contrast to the results obtained from polarizable potentials, the corresponding nonpolarizable number density profile had a depletion of both Cl- and Sr2+ ions, with all the ions concentrated in the bulk of solution (results not shown). This is a similar observation to what was found from simulations of monovalent ions with nonpolarizable potential.11,19,39,40 The total electron density profiles obtained from simulations with both nonpolarizable and polarizable potentials are shown in Figure 2 (see inset). The curves were calculated from an average of the two liquid/vapor interfaces. The center of the bulk solution is located at 0 Å in this figure. The electron density profiles from simulations with nonpolarizable and polarizable potentials show tremendous differences. The calculated bulk electron density is 0.37 e/Å3 for a uniformly distributed ion solution (shown as blue dashed line). The computed electron density from the simulations with polarizable potentials has a peak height of 0.53 e/Å3 at the surface and drops to a bulk value of 0.35 ( 0.01 e/Å3 about 7 Å away from the surface inside the solution. The high electron density observed at the surface is clearly from the high concentration of ions at the surface, and especially Sr2+ ions, which have a much higher atomic number but smaller molecular size due to the fact that it is a divalent cation, giving it a very high electron density compared to water molecules and Cl-. The electron density computed from the simulations with nonpolarizable potentials shows a monotonic decrease from liquid to the vapor phase with a bulk value of 0.41 ( 0.01 e/Å3. The contribution to this significant difference between electron densities is that the ions are depleted from the surface for the nonpolarizable models and concentrated at the surface for the polarizable models.
Aqueous Interfaces of SrCl2 Salt Solutions
J. Phys. Chem. B, Vol. 113, No. 42, 2009 13995
Figure 2. Comparison of computed |Φ(qz)|2 from nonpolarizable (red) and polarizable (green) potential models and experimentally measured values (symbols). (Inset) Comparison of simulated total electron density profiles for nonpolarizable and polarizable potential models. The blue dashed line is the electron density profile for uniformly distributed ions in aqueous solution.
As mentioned in our previous work, the structure factor is a suitable property to depict the nature of the surface structure.19 Figure 2 shows the calculated |Φ(qz)|2 from simulations with nonpolarizable and polarizable potentials as well as the experimental results. Because the simulation environment is substantially different from the experimental one, such as the system size, equilibration time, etc., the calculated |Φ(qz)|2 values from the simulations are optimized by a Debye-Waller factor to scale the amplitude to be a better fit to the experiments. Generally, this extra term is small and will not change the resonance frequency of the structure factor in our calculation. The experimentally calculated |Φ(qz)|2 from X-ray reflectivity measurements shows a single peak at qz ≈ 0.4 Å-1 and large fluctuations in the high-frequency region. The calculated |Φ(qz)|2 from the simulations with polarizable potentials qualitatively agree with the experimental curve. The same resonance frequency is observed at qz ≈ 0.4 Å-1, and the trend of the |Φ(qz)|2 corresponding to the frequency follows the experiment reasonably well. However, the peak of simulated |Φ(qz)|2 is more abrupt than that of experiment. The simulations with nonpolarizable potentials show a much different scenario. Although an experimentally matched resonance frequency can be observed, the shape of the curve is significantly different from the experimental one. This could not be further improved by any choice of Debye-Waller factor. The comparison between the |Φ(qz)|2 calculated from simulations and X-ray reflectivity indicates an enhanced distribution of ions at the liquid/vapor interface. However, the simulations with polarizable potentials give overstructured spectra. In ref 27, since there are no direct MD simulations of SrCl2 solutions at the experimental concentration, an electron density profile generated from a lower concentrated NaCl solution simulation was used to calculate the structure factor, and no agreement with experiment was found. Our simulations with polarizable potentials provide a direct comparison with experiment for the interfacial structural properties at the molecular level. It is possible that a lower interfacial ion concentration than what is found from our simulations with polarizable potentials would improve agreement with experiment. Further efforts on improving the ion-ion and ion-water interactions should be addressed in the future. To better understand the effect of the electron density profile on the structure factor, two tests were performed. Since Sr2+
Figure 3. (a) Results of the computed |Φ(qz)|2 when the contributions of Sr2+ to the total electron density were artificial reduced to 0%, 50%, and 100%. (b) Simulated |Φ(qz)|2 as a function of SrCl2 salt concentrations. (c) Same as Figure 1 but for 1.35 M SrCl2.
has a relatively small ionic radius and large atomic number, the electron density of Sr2+ is much higher than either a Clatom or a water molecule, and it is expected that the existence of Sr2+ has a large impact on the total electron density profiles. If the contribution of Sr2+ to the total electron density is artificially reduced by 50% and 100%, the corresponding ratios of surface to bulk electron density (rsurface/rbulk) will be reduced from 1.5 to 1.2 and 1.0, respectively. The calculated |Φ(qz)|2 with these assumptions are shown in Figure 3a. Considerable changes can be observed. The resonance peak at qz ≈ 0.4 Å-1 is lowered when the contribution of Sr2+ to the total electron density profile is reduced. When the contribution is reduced by 50%, the calculated |Φ(qz)|2 agrees with the experimental curve very well. When the Sr2+ contribution is reduced by 100%, the total electron density profile is almost monotonic and the resonance peak disappears. The tests indicate that Sr2+ has a non-negligible contribution to the total electron density profile at the surface. The resonance peak at qz ≈ 0.4 Å-1 appears to result from enhancement of the total electron density at the liquid/vapor interface, and the amplitude of the peak is dominated by the degree of this enhancement. Although the agreement between calculated |Φ(qz)|2 from the 50% Sr2+
13996
J. Phys. Chem. B, Vol. 113, No. 42, 2009
Sun et al. that multiple sets of these parameters can be used to give similar |Φ(qz)|2 values in the low-frequency region. Our simulation results provide a direct comparison to be made with X-ray reflectivity experiments and allow the spectroscopic features to be linked to interfacial structure and ion distributions. IV. Conclusion
Figure 4. Simulated |Φ(qz)|2 from simulations with polarizable SrCl2 (green) and NaCl (red) potential models.
contribution and that of experiment is only a numerical coincidence, this test provides some insights into the experimental electron density profile. The experimental ratio of surface to bulk electron density is suggested to be about 1.2. To understand the role of salt concentration, we carried out additional simulations at lower concentrations (1.00 and 1.35 M) of SrCl2 salt solutions. The simulated X-ray reflectivity results are presented in Figure 3b. Upon examination of these results, better agreement with experiment can be achieved with the 1.35 M SrCl2 (50% of 2.7 M) salt solution. It is also very similar to the result in which the Sr2+ contribution was artificially reduced by 50%. In Figure 3c, we present the calculated number density and a snapshot of the simulation of 1.35 M SrCl2. We find the results are very similar to the results of the 2.7 M SrCl2 aqueous salt solution shown in Figure 1. We can conclude that the origin of the overstructured X-ray spectra calculated from our MD simulations is likely the result of the surface concentration of SrCl2 ions being around a factor of 2 too high. This was confirmed both by the improved agreement between simulation and experiment when the Sr2+ contribution to the spectra was halved in a 2.7 SrCl2 solution and by the agreement shown when the overall SrCl2 concentration is halved. However, the results also show that there appear to be a significant number of Sr2+ ions present at the surface, which is in contrast to what is observed for the monovalent cations (which form a double layer). Figure 4 shows |Φ(qz)|2 calculated from simulations of a 2.7 M NaCl solution. Due to the lower valency of Na+ compared to that of Sr2+, the ions spread in the solution as solvent-separated ion pairs. A weak enhancement of Cl- was observed at the liquid/vapor interface, and Na+ is enhanced several angstroms away from the interface where Cl- shows a decrease (not shown). The ratio of the surface to bulk electron density is about 1.1. Although the total electron density distribution is very different from that of SrCl2, a similar quality fit to experimentally calculated |Φ(qz)|2 of SrCl2 solution as the fit from simulated 50% Sr2+ contributed SrCl2 solution can be observed in the lower and midrange frequency region. However, unlike the 50% Sr2+ contributed SrCl2 solution and the solutions with lower SrCl2 concentrations, deviation of calculated |Φ(qz)|2 from NaCl with respect to that of experiment becomes more divergent at higher frequency. This shows that solutions with NaCl and SrCl2 share some similar features but give fundamentally different spectra. A single-slab model has been implemented to investigate the relationship between |Φ(qz)|2 and electron density profile.27 In this model, three free parameterssd, the thickness of the surface layer; Fs, the surface electron density; and σ0, the intrinsic roughnessshave being used to characterize the liquid/vapor interface but have been shown to be interdependent. This means
SrCl2 aqueous solutions were studied by molecular dynamics simulations with both nonpolarizable and polarizable potential models and at various concentrations. The results show that both Sr2+ and Cl- are enhanced at the liquid/vapor interface from the simulations with polarizable potential models, and both are depleted from the interface from simulations with nonpolarizable potential models. The structure factor calculated from the simulations with nonpolarizable potential models cannot reproduce the experimentally calculated one, even qualitatively. The simulations with polarizable potential models give qualitative agreement with experimental measured data. Furthermore, if the contribution of Sr2+ to the structure factor is reduced by 50%, near-quantitative agreement with experiment can be obtained. We also performed simulations at lower salt concentrations (1.0 and 1.35) using polarizable potential models, and significant improvement with experimental spectra can be achieved when at 1.35 M. This led us to conclude that the surface concentration of SrCl2 is significant, but our molecular models overestimate the surface concentration by around a factor of 2. In addition, this study demonstrated that simulations of NaCl aqueous solution gave significantly different results with the corresponding simulation results of SrCl2 X-ray reflectivity. We can attribute this finding to the fact that, unlike the monovalent ions, which form solvent-separated ion pairs in dilute solutions, Sr2+ and Cl- primarily exist as contact ion pairs in the solution due to the large net charge on Sr2+. Acknowledgment. This work was performed at Pacific Northwest National Laboratory under the auspices of the Division of Chemical Sciences, Geosciences, and Biosciences, Office of Basic Energy Sciences, U.S. Department of Energy (DOE). Pacific Northwest National Laboratory is operated by Battelle for the DOE. The Division of Chemical and Materials Sciences provided computer resources. References and Notes (1) Jones, G.; Ray, W. A. J. Am. Chem. Soc. 1937, 59, 187. (2) Dole, M.; Swartout, J. A. J. Am. Chem. Soc. 1940, 62, 3039. (3) Weissenborn, P. K.; Pugh, R. J. J. Colloid Interface Sci. 1996, 184, 550. (4) Markin, V. S.; Volkov, A. G. J. Phys. Chem. B 2002, 106, 11810. (5) Gibbs, J. W. The Collected Works of J. Willard Gibbs; Longmans: New York, 1928. (6) Langmuir, I. J. Am. Chem. Soc. 1917, 39, 1848. (7) Wagner, C. Phys. Z. 1924, 25, 474. (8) Onsager, L.; Samaras, N. N. T. J. Chem. Phys. 1934, 2, 528. (9) Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2001, 105, 10468. (10) Dang, L. X. J. Phys. Chem. B 2002, 106, 10388. (11) Dang, L. X.; Chang, T. M. J. Phys. Chem. B 2002, 106, 235. (12) Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2002, 106, 6361. (13) Garrett, B. C. Science 2004, 303, 1146. (14) 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. (15) Chang, T. M.; Dang, L. X. Chem. ReV. 2006, 106, 1305. (16) Jungwirth, P.; Tobias, D. J. Chem. ReV. 2006, 106, 1259. (17) Wick, C. D.; Dang, L. X. J. Phys. Chem. B 2006, 110, 6824. (18) Wick, C. D.; Dang, L. X. Chem. Phys. Lett. 2008, 458, 1. (19) Sun, X. Q.; Dang, L. X. J. Chem. Phys. 2009, 130, 124709. (20) 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. (21) Liu, D. F.; Ma, G.; Levering, L. M.; Allen, H. C. J. Phys. Chem. B 2004, 108, 2252.
Aqueous Interfaces of SrCl2 Salt Solutions (22) Petersen, P. B.; Johnson, J. C.; Knutsen, K. P.; Saykally, R. J. Chem. Phys. Lett. 2004, 397, 46. (23) Petersen, P. B.; Saykally, R. J. Chem. Phys. Lett. 2004, 397, 51. (24) Ghosal, S.; Hemminger, J. C.; Bluhm, H.; Mun, B. S.; Hebenstreit, E. L. D.; Ketteler, G.; Ogletree, D. F.; Requejo, F. G.; Salmeron, M. Science 2005, 307, 563. (25) Petersen, P. B.; Saykally, R. J. J. Phys. Chem. B 2005, 109, 7976. (26) Petersen, P. B.; Saykally, R. J.; Mucha, M.; Jungwirth, P. J. Phys. Chem. B 2005, 109, 10915. (27) Sloutskin, E.; Baumert, J.; Ocko, B. M.; Kuzmenko, I.; Checco, A.; Tamam, L.; Ofer, E.; Gog, T.; Gang, O.; Deutsch, M. J. Chem. Phys. 2007, 126, 054704. (28) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press Inc.: New York, 1987. (29) Brown, E. C.; Mucha, M.; Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2005, 109, 7934. (30) Ishiyama, T.; Morita, A. J. Phys. Chem. C 2007, 111, 738.
J. Phys. Chem. B, Vol. 113, No. 42, 2009 13997 (31) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (32) Dang, L. X.; Chang, T. M. J. Chem. Phys. 1997, 106, 8149. (33) Smith, D. E.; Dang, L. X. Chem. Phys. Lett. 1994, 230, 209. (34) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757. (35) Ulrich, E.; Lalith, P.; Max, L. B.; Tom, D.; Hsing, L.; Lee, G. P. J. Chem. Phys. 1995, 103, 8577. (36) Berendsen, H. J. C.; Postma, J. P. M.; Gunsteren, W. F. v.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (37) Sloutskin, E.; Lynden-Bell, R. M.; Balasubramanian, S.; Deutsch, M. J. Chem. Phys. 2006, 125, 174715. (38) Friedman, H. L. Ionic Solution Theory, Based on Cluster Expansion Methods; Interscience Publishers: New York, 1962. (39) Benjamin, I. J. Chem. Phys. 1991, 95, 3698. (40) Wilson, M. A.; Pohorille, A. J. Chem. Phys. 1991, 95, 6005.
JP9079525
