Journal of Molecular Liquids 112 (2004) 47–50
Charge separation at the iceywater interface: a molecular dynamics simulation study of solute ions at the ice basal plane Taras Bryka,c, A.D.J. Haymetb,c,* a
Institute for Condensed Matter Physics, National Academy of Sciences of Ukraine, 1 Svientsitskii Street, UA-79011 Lviv, Ukraine b CSIRO Marine Research, GPO Box 1538, Hobart TAS 7001 Australia c Department of Chemistry, University of Houston, Houston, TX 77204-5003 USA
Abstract The behavior of Naq and Cly ions over a total time of 1.5 ns ice basal planeyvacuum interface is simulated with the SPCyE model of water molecules at a temperature close to the coexistence point. It is found that the Naq ion stays for almost the whole simulation time in the top broad ice surface layer, while Cly ion penetrates through the top layer and settles in the subsurface layer. Running coordination numbers calculated every 60 ps show on average seven nearest neighbors for the Cly ion, while the Naq ion is surrounded mostly by five nearest water molecules. A tendency to melt is observed for the surface ice layers with the solute ions inside. 䊚 2003 Elsevier B.V. All rights reserved. PACS: 82.65.qr; 68.35.-p; 82.20.Wt Keywords: Charge separation; Ice surface; Ice–water interface; Computer simulations
1. Introduction Studies of the iceywater interface are of great interest because of their fundamental presence in nature and importance in chemical, biological, environmental and atmospheric processes w1x. Many fundamental phenomena occurring on the ice surface and iceywater interface still have not been explained on a microscopic level. More than 50 years ago, Workman and Reynolds w2x claimed to observe a massive charge separation occurring during the freezing of dilute aqueous solutions, and imputed a relevance to thunderstorm electricity. These experiments have never been duplicated successfully. The observed transient and then static effect has not been explained to date, notwithstanding an attempt to explain the kinetic effect using the different probabilities of positive and negative ions trapped into ice during the freezing front propagation w3x. To the best of our knowledge, there are also no computer simulation studies with realistic models, which seek microscopic insight into processes occurring with ions close to the iceywater interface. Additional interest for this study arises from *Corresponding author. E-mail address:
[email protected] (A.D.J. Haymet).
alleged catalytic properties of the interfaces. It is supposed that the HCl and HBr acids can dissociate at the ice surface w4–6x due to the presence of top liquid-like layers on the ice surface. An essential advantage of computer simulations over phenomenological approaches is the possibility of observing the microscopic processes on atomic time and spatial scales. In this study, we investigate the behavior of Naq and Cly ions on the ice surface, which is modeled using a realistic model of water molecules. In molecular dynamics (MD) simulations, the solute ions move from an initial position towards a final one, which is characterized by a smaller free energy of the whole system. Thus, we expect in this study to observe the solute ion evolving towards an equilibrium state with lower free energy. 2. Simulation method We have simulated the basal ice surface with Naq and Cly ions in the NVT ensemble at the temperature of 225 K, which was estimated as the melting point for SPCyE model w7x of water in our previous study of icey water coexistence w8x. The basal surface of ice was prepared from a well-equilibrated ice Ih bulk phase
0167-7322/04/$ - see front matter 䊚 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0167-7322(03)00274-5
48
T. Bryk, A.D.J. Haymet / Journal of Molecular Liquids 112 (2004) 47–50
3. Results and discussion 3.1. Analysis of ion location on the ice basal surface
Fig. 1. Mass-density (upper frame) and charge-density profiles for the equilibrium basal ice surface at 225 K without solute ions.
simulated with 2304 SPCyE molecules. A middle part of the ice crystal was cut-off creating a vacuum region ˚ We allowed the in MD box with the width of ;37 A. ice surface to relax, and then observed it over the time of ;300 ps in order to make sure that it was stable. The mass-density and charge-density profiles of the ice basal surface are shown in Fig. 1. One can observe the specific double peaks in the bulk ice Ih phase on the ˚ while the presence of mass-density profile for z)29 A, the surface is reflected in smearing the peaks for z-29 ˚ The top layer does not have the double peak structure A. anymore due to strong smearing. The charge-density profile permits the analysis of the regions with excess of hydrogens or oxygens. One can see that the water ˚ molecules on the top of the ice surface, 17 AFzF19 ˚ A, have the preferential orientation with oxygens headed towards the bulk ice phase. The charge-density profile itself implies the possibility for solute ions with opposite charges to be localized in different depths. In these simulations the time step was 1.5 fs, the same as in the previous study w8x. The details on shortrange interaction and Ewald summation can be found also elsewhere w8x. For Lennard–Jones interaction between the ions and SPCyE water molecules, we have chosen the parameters proposed originally by Smith and Dang w9x and used later in simulations by Kovalenko and Hirata w10x. All MD simulations were performed by the standard DL_POLY package w11x.
˚ in the MD The ions were initially placed at zs15 A ˚ from the box, i.e. approximately at the distance of 2 A surface, then we allowed the ions to approach the interface during the time of 20 ps, and only after that started the production runs. In order to characterize the location of solute ions with respect to the surface we introduced a single-ion density profile accumulated over a certain time window. We have found that such singleion density profiles for the case of our system, being accumulated over the time of 60 ps, have the width of ˚ and contain a well-defined peak, which corref2 A sponds to the preferential location of the ion with respect to the surface over the 60 ps. In Fig. 2, we show examples of single-ion density profiles for Naq (upper frame) and Cly (lower frame). One can see that some of the single-density profiles can have shoulders and small local maxima. They correspond to the local minima of potential energy surface, at which the ions spend several picoseconds oscillating around the minimum location. These oscillations can be easily observed on a z-component of force autocorrelation functions, introduced in our recent simulations of ions at the basal iceywater interface w12x. The change in the main maximum location reflects the location of the ions with respect to the ice surface.
Fig. 2. Single-ion density profiles for Naq (upper frame) and Cly ions, accumulated during different time windows.
T. Bryk, A.D.J. Haymet / Journal of Molecular Liquids 112 (2004) 47–50
Fig. 3. Location of Naq ion on the ice basal surface over the simulation time of 1.5 ns (symbols). Each symbol ‘q’ corresponds to the position of maximum in a single-ion density profile accumulated over each time window of 60 ps. Solid line in the right frame corresponds to the mass-density profile accumulated over the first 60 ps.
In Figs. 3 and 4 one can see that the Naq and Cly ions behave in different way being on the ice surface. Naq is moving within the top strongly smeared double layer, ˚ ˚ while the Cly ion penetrates through 20 AFzF22 A, the top surface layer during the first 120 ps and settles on the top side in the second to the surface double layer, ˚ ˚ (see Fig. 4). Figs. 3 and 4 are evidence 22 AFzF24 A that the ice surface simulated with the SPCyE model of rigid water molecules can produce a separation of oppositely charged ions with respect to the surface. To find an origin of the charge separation in the case of Naq and Cly ions on the basal ice surface, we will analyze the running coordination numbers nNa–O(r) and nCl–O(r), obtained from the relevant radial distribution functions during each time window. 3.2. Running coordination numbers Since the center of mass for water molecules is located very close to the oxygen position, the analysis
Fig. 4. Location of Cly ion on the ice basal surface over the simulation time of 1.5 ns (symbols). Each symbol ‘q’ corresponds to the position of maximum in a single-ion density profile accumulated over each time window of 60 ps. Solid line in the right frame corresponds to the mass-density profile accumulated over the first 60 ps.
49
Fig. 5. Running coordination numbers for oxygens around Naq ion on the ice basal surface during different time windows. The time window number corresponds to the labels in Fig. 3.
of running coordination numbers ion–oxygen reflects the molecular environment around the ion. In Fig. 5, one can see that the Naq ion is surrounded by five nearest water molecules on the ice basal surface. One should note that the five-molecule first hydration shell is possible due to the defects in the SPCyE ice structure caused by the Naq ion. For constrained (rigid) SPCyE ice lattice the Naq ion would have a local minimum in a plane of hexagonal hole of the ice structure with the six nearest water molecules. Although in Fig. 5 only the time windows with five-molecule first hydration shell are shown, we also observed the possibility of sixmolecule hydration shells for some other time windows. We observed a tendency to melting the surface ice structure with the Naq ion inside. In Fig. 6, the massdensity profiles are compared at three different time windows. Obvious tendency to a liquid-like region is reflected in a very strong smearing of the top and next to the top ice double layers. We observed an increase in the time spent by the Naq ion in the six-molecule hydration shells during such a melting of the ice surface.
Fig. 6. Mass-density profile of the ice basal surface with the Naq ion during different time windows. The time window number corresponds to the labels in Fig. 3.
50
T. Bryk, A.D.J. Haymet / Journal of Molecular Liquids 112 (2004) 47–50
4. Conclusions
Fig. 7. Running coordination numbers for oxygens around Cly ion on the ice basal surface during different time windows. The time window number corresponds to the labels in Fig. 4.
The Cly ions cause much stronger defects in the ice structure of SPCyE water molecules. Repulsion between oxygens and Cly ion leads to an average radius of the ˚ while for the confirst hydration shell of ;3.25 A, strained (rigid) SPCyE ice lattice it would be smaller, ˚ Also, the Cly ions have a larger core, which ;2.82 A. is reflected in higher coordination numbers for the ion on the ice surface. In Fig. 7, one can see that the Cly ion has seven SPCyE water molecules in the first hydration shell. But more interesting is a change in the second hydration shell for the Cly ion. Increasing the molecules in the second hydration shell is possible for the ion only when it settles far from the surface. This can be a reason why the Cly ion is located more deep under the ice basal surface than the Naq ion. In Fig. 8, the similar tendency to melting the surface structure with Cly inside is shown. We can conclude that both ions, Naq and Cly, favor the ice basal surface melting at the temperatures close to the melting point in the model of SPCyE water molecules.
In this study, we investigated by MD computer simulations the behavior of Naq and Cly ions at the ice basal surface using the SPCyE rigid model of water molecules. We have found that the ions, initially placed ˚ from the surface, evolve over a at the distance of 2 A total time of approximately 1.5 ns at different depths inside the interface. The Naq ion stays in the top surface layer, while the Cly ion penetrates the top layer during first 120 ps and then settles in the subsurface layer. Analysis of the running coordination numbers implies that the different equilibrium location of the ions inside the ice basal surface is defined by a higher number of neighbors for Cly, which may be the consequence of repulsion between oxygens and Cly ion, and a size effect, namely the larger radius of the Cly ion compared to Naq. Another interesting finding is the tendency to melt for the surface and subsurface ice layers containing the solute ions. Since we have used only the SPCyE model of rigid water molecules in our simulations, it would be reasonable to undertake similar simulations based on a more realistic model of polarizable water molecules in order to check both effects found in this study, i.e. the charge separation and the tendency to melting of top layers with solute ions inside. Acknowledgments This research was supported by the Welch Foundation Grant E-1429. The calculations were performed on CrayT3E and IBM-Power4 at Texas Advanced Computing Center under NPACI award. DL_POLY is a package of molecular simulation routines written by W. Smith and T.R. Forester, copyright The Council for the Central Laboratory of the Research Councils, Daresbury Laboratory at Daresbury, Nr. Warrington (1996). References
Fig. 8. Mass-density profile of the ice basal surface with the Cly ion during different time windows. The time window number corresponds to the labels in Fig. 4.
w1x V.F. Petrenko, R.W. Whitworth, Physics of Ice, University Press, Oxford, 1999. w2x E.J. Workman, S.E. Reynolds, Phys. Rev. 78 (1950) 254. w3x V.L. Bronshteyn, A.A. Chernov, J. Cryst. Growth. 112 (1991) 129. w4x B.J. Gertner, J.T. Hynes, Science 271 (1996) 1563. w5x D.C. Clary, L. Wang, Faraday Trans. 93 (1997) 2763. w6x M. Svanberg, J.B. Pettersson, K. Bolton, J. Phys. Chem. A 104 (2000) 5787. w7x H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. w8x T. Bryk, A.D.J. Haymet, J. Chem. Phys. 117 (2002) 10 258. w9x D.E. Smith, L.X. Dang, J. Chem. Phys. 100 (1994) 3757. w10x A. Kovalenko, F. Hirata, J. Chem. Phys. 112 (2000) 10 391. w11x http:yywww.dl.ac.ukyTCSySoftwareyDL_POLYy. w12x E.J. Smith, T. Bryk, A.D.J. Haymet, (in preparation).