Molecular Dynamics Simulation of Aqueous Sodium Chloride Solution

More than 80% of all the potassium minerals in the world are recovered in this way (as of the early 1990s).16 Moreover, interfaces between sodium chlo...
22 downloads 0 Views 145KB Size
Download PDF
Langmuir 2002, 18, 547-556

547

Molecular Dynamics Simulation of Aqueous Sodium Chloride Solution at the NaCl(001) Interface with a Polarizable Water Model Enno Oyen and Reinhard Hentschke* FB Physik, Bergische Universita¨ t-Gesamthochschule Wuppertal, D-42097 Wuppertal, Germany Received August 10, 2001. In Final Form: October 25, 2001 The sodium chloride solution-NaCl crystal interface serves as an example of aqueous electrolytes at uncharged ionic surfaces and plays an important role in technologically relevant processes and atmospheric chemistry. This interface is investigated by classical molecular dynamics simulations at ambient pressure and generally at 298 K. A fluctuating charges version of the extended simple point charge (SPC/E) water model, the SPC/E-P model, is used to calculate structural and dynamical properties of the ions at the interface and in the bulk solution: density profiles and potentials of mean force including concentration and temperature dependencies, hydration numbers, lateral structures, and mobilities of ions parallel and perpendicular to the interface. We find that sodium ions are adsorbed selectively at the metastable NaCl(001) interface, and their mobility is strongly reduced due to adsorption. Analogous calculations using the nonpolarizable SPC/E model do not exhibit such an intense exothermic sodium ion adsorption; i.e., solvent polarization can significantly affect ionic distributions close to polar but uncharged interfaces. In the case of SPC/E-P water, dipole moment profiles are calculated also.

1. Introduction Modeling electrical double layers at interfaces between metals and aqueous electrolytes has been and still is an active area of research (for a recent review see ref 1). Investigating water or aqueous solutions at mineral surfaces is less common, but a few examples of molecular dynamics (MD) simulations using phenomenological potentials include aqueous phases at magnesium oxide,2 the aluminosilicate kaolinite,3,4 and sodium chloride. In the latter case, classical MD simulations of the saltwater interface range from water adsorption on NaCl crystals5,6,7 to bulk aqueous phases in contact with NaCl crystals.7-10 Quantum-mechanical treatments of NaCl(001) interfaces with water comprise semiempirical, HartreeFock level, and embedded cluster density functional calculations of water adsorption on clean surfaces11-13 as well as an embedded cluster density functional approach to the bulk water-salt interface with adsorption of water molecules and sodium and chloride ions.14 * To whom correspondence should be addressed. (1) Guidelli, R.; Schmickler, W. Electrochim. Acta 2000, 45, 23172338. (2) Marmier, A.; Hoang, P. N. M.; Picaud, S.; Girardet, C.; LyndenBell, R. M. J. Chem. Phys. 1998, 109, 3245-3254. (3) Mattke, T.; Kecke, H.-J. J. Colloid Interface Sci. 1998, 208, 555561. (4) Mattke, T.; Kecke, H.-J. J. Colloid Interface Sci. 1998, 208, 562569. (5) Wassermann, B.; Mirbt, S.; Reif, J.; Zink, J. C.; Matthias, E. J. Chem. Phys. 1993, 98, 10049-10060. (6) Sto¨ckelmann, E.; Aydt, E. M.; Hentschke, R. J. Mol. Model. 1997, 3, 347-354. (7) Sto¨ckelmann, E.; Hentschke, R. J. Chem. Phys. 1999, 110, 1209712107. (8) Ohtaki, H.; Fukushima, N. Pure Appl. Chem. 1989, 61, 179-185. (9) Shinto, H.; Sakakibara, T.; Higashitani, K. J. Phys. Chem. B 1998, 102, 1974-1981. (10) Shinto, H.; Sakakibara, T.; Higashitani, K. J. Chem. Eng. Jpn. 1998, 31, 771-779. (11) Jug, K.; Geudtner, G. Surf. Sci. 1997, 371, 95-99. (12) Taylor, D. P.; Hess, W. P.; McCarthy, M. I. J. Phys. Chem. B 1997, 101, 7455-7463. (13) Stefanovich, E. V.; Truong, T. N. J. Chem. Phys. 1996, 104, 29462955.

Aqueous sodium chloride solutions at NaCl crystal surfaces serve as examples for aqueous electrolytes in contact with ionic interfaces. Furthermore, they are relevant on their own. The surface charges (especially their signs) of alkali halide particles in their saturated solutions, which result from selective ion adsorption, play an important role in processes such as corrosion, tertiary oil recovery, crystallization, and flotation.15 Flotation is used on an industrial scale for separating soluble salts such as NaCl and KCl: The suspension of hydrophilic salt particles in their saturated solution is mixed with a collector, a surfactant selectively binding to the particles of a single salt and thus turning their surfaces into a hydrophobic state. Air bubbles dispersed in the suspension raise the hydrophobic particles attached to them to the surface, where these particles can be removed separately. Miller et al.15 showed that the crucial adsorption of charged collectors can be explained in terms of Coulombic interactions between collector charges and surface charges of alkali halide particles. More than 80% of all the potassium minerals in the world are recovered in this way (as of the early 1990s).16 Moreover, interfaces between sodium chloride and water or aqueous solutions play an important role in environmental chemistry, because atmospheric sea salt particles with NaCl as major constituent sometimes form crystals covered with water rather than concentrated aqueous solutions or anhydrous minerals.17 The water layer can alter reactions of NaCl surfaces with pollutants such as nitrous gases and nitric acid.18 Sea salt particles are probably the dominant source of reactive atmospheric chlorine compounds, which affect the sensitive ozone balance and oxidize organic substances.19 (14) Stefanovich, E. V.; Truong, T. N. J. Chem. Phys. 1997, 106, 77007705. (15) Miller, J. D.; Yalamanchili, M. R.; Kellar, J. J. Langmuir 1992, 8, 1464-1469. (16) Veeramasuneni, S.; Yalamanchili, M. R.; Miller, J. D. J. Colloid Interface Sci. 1996, 182, 275-281. (17) Allen, H. C.; Laux, J. M.; Vogt, R.; Finlayson-Pitts, B. J.; Hemminger, J. C. J. Phys. Chem. 1996, 100, 6371-6375. (18) Vogt, R.; Finlayson-Pitts, B. J. J. Phys. Chem. 1995, 99, 1726917272.

10.1021/la011269e CCC: $22.00 © 2002 American Chemical Society Published on Web 12/27/2001

548

Langmuir, Vol. 18, No. 2, 2002

There are numerous papers in the literature dealing with modeling homogeneous sodium chloride solutions. The methods used comprise Monte Carlo simulations,20,21 molecular dynamics,22-28 and calculations based on integral equation theory.29-31 Rigid water models without polarization such as the four-point transferable intermolecular potential (TIP4P),32 the simple point charge (SPC) model,33 and the extended simple point charge (SPC/E) model34 are employed in most of these works. Only a few authors make use of flexible or polarizable models. Although homogeneous systems are simpler to simulate than inhomogeneous solutions at interfaces, there are long-standing, yet unresolved, problems, e.g., the ion-ion radial distribution functions (RDFs) in aqueous NaCl solutions. These RDFs cannot be obtained experimentally since neutron scattering data are not sufficiently accurate.31 Furthermore, experimental data concerning the structure of aqueous sodium ions is limited to the first hydration shell.35 So molecular simulation is the only available method for obtaining ion-ion RDFs in electrolyte solutions at moderate concentrations. There are technical and conceptual problems, however: On the one hand, quite long simulation times of several nanoseconds are required for calculating ion-ion RDFs from molecular dynamics simulations of aqueous NaCl solutions due to the small number and slow diffusion of ions.25,26 On the other hand and even worse, the resulting ion-ion RDFs of NaCl solutions depend on the detail of potential models employed for water and ions.26,31 An example of a molecular dynamics study of an inhomogeneous system is the one by Kohlmeyer et al.,36 who modeled a water/metal interface using both a polarizable and a nonpolarizable water model. They found that atomic density profiles of their models barely depend on polarization and that the dipole moment of polarizable water at the wall is only slightly decreased as compared to the bulk phase, though the hydration layers are much more dense than bulk water. The authors arrived at the conclusion that polarization is not needed for modeling metal-water interfaces. This conclusion is supported by (19) Laux, J. M.; Fister, T. F.; Finlayson-Pitts, B. J.; Hemminger, J. C. J. Phys. Chem. 1996, 100, 19891-19897. (20) Degre`ve, L.; da Silva, F. L. B. J. Chem. Phys. 1999, 110, 30703078. (21) Lyubartsev, A. P.; Laaksonen, A. Phys. Rev. E 1997, 55, 56895696. (22) Degre`ve, L.; da Silva, F. L. B. J. Chem. Phys. 1999, 111, 51505156. (23) Driesner, T.; Seward, T. M.; Tironi, I. G. Geochim. Cosmochim. Acta 1998, 62, 3095-3107. (24) Hermansson, K.; Wojcik, M. J. Phys. Chem. B 1998, 102, 60896097. (25) Koneshan, S.; Rasaiah, J. C. J. Chem. Phys. 2000, 113, 81258137. (26) Lyubartsev, A. P.; Laaksonen, A. J. Phys. Chem. 1996, 100, 16410-16418. (27) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757-3766. (28) Zhu, S.-B.; Robinson, G. W. J. Chem. Phys. 1992, 97, 43364348. (29) Fedotova, M. V.; Trostin, V. N. Russ. J. Phys. Chem. 1999, 73, 905-908. (30) Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 1039110402. (31) Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 1040310417. (32) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (33) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces: Proceedings of the Fourteenth Jerusalem Symposium on Quantum Chemistry and Biochemistry; Pullman, B., Ed.; Reidel: Dordrecht, 1981; pp 331-342. (34) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269-6271. (35) Kim, H.-S. Chem. Phys. 2000, 253, 305-312. (36) Kohlmeyer, A.; Witschel, W.; Spohr, E. Chem. Phys. 1996, 213, 211-216.

Oyen and Hentschke

comparing results from MD simulations of water near NaCl(001) interfaces employing nonpolarizable SPC/E water9 and the polarizable SPC/E-P model,7 respectively: Several hydration layers of immobilized water molecules and similar atomic density profiles are observed in both cases, and the dipole moment of polarizable water at the interface is only slightly smaller than its bulk value. Nevertheless, in MD simulations of nonpolarizable water at an NaCl(001) surface, the adsorption of water turned out to be strongly affected by the exact value of the dipole moment.6 Hence, Sto¨ckelmann et al. argue for explicitly including polarization when modeling water in strongly inhomogeneous electric fields close to ionic interfaces. Moreover, Stuart et al.37 observed a clear effect of water polarization, but no effect of chloride polarization on the solvation of chloride ions in model water clusters. The water structure at a wall can also influence or even determine the distribution of ions.38,39 Even differences between density profiles of sodium and chloride ions, which were found in simulations of aqueous electrolyte at a metal electrode with identical metal-ion potentials for cations and anions, could be attributed to differences in ionic hydration alone.40 Comparing results obtained with two simple hard sphere electrolyte models, the restricted primitive and the molecular solvent model, at a wall, Patra et al.41-43 showed that layering is intensified and ionic density profiles oscillate considerably if solvent particles are included explicitly instead of using a continuum solvent model, which also emphasizes the importance of solvent effects. Thus, currently one may argue both in favor of and against explicitly considering solvent polarization. Consequently, we decided to employ an efficient polarizable water model in our MD simulations of the aqueous salt solution-NaCl(001) interface and to directly compare the results with those obtained from a similar nonpolarizable model. Hydration layers have already been observed in oxygen density profiles from X-ray scattering experiments of charged Ag(111)-electrolyte interfaces.44 On the contrary, there is still little known from experiment about ionic distributions in electrolytes at walls. Fenter et al.,45 who carried out surface X-ray diffraction and X-ray adsorption spectroscopy at the rutile(110)-electrolyte interface with adsorbed rubidium and strontium ions, managed to measure ion positions within the condensed part of the electrical double layer and the ion partitioning between the condensed and diffuse layers, but no ion density profile. Ion density profiles can, therefore, be obtained from simulation and other methods of modeling only. As early as 1989, Ohtaki et al.8 simulated a sodium chloride crystal dissolving in water leading to a positively charged crystal, because several chloride, but no sodium, ions were detached from the crystal within the extremely short period of 7 ps covered by their MD run. Analyzing this kinetic effect, the authors traced back the expulsion of chloride ions from the crystal to the repulsion of chloride by water molecules attracted by crystal sodium ions. We (37) Stuart, S. J.; Berne, B. J. J. Phys. Chem. 1996, 100, 1193411943. (38) Be´rard, D. R.; Kinoshita, M.; Cann, N. M.; Patey, G. N. J. Chem. Phys. 1997, 107, 4719-4728. (39) Vossen, M.; Forstmann, F. Mol. Phys. 1995, 86, 1493-1516. (40) Spohr, E. Electrochim. Acta 1999, 44, 1697-1705. (41) Patra, C. N. J. Chem. Phys. 1999, 111, 9832-9838. (42) Patra, C. N.; Ghosh, S. K. J. Chem. Phys. 1994, 101, 4143-4149. (43) Patra, C. N.; Ghosh, S. K. J. Chem. Phys. 1994, 100, 5219-5229. (44) Toney, M. F.; Howard, J. N.; Richer, J.; Borges, G. L.; Gordon, J. G.; Melroy, O. R.; Wiesler, D. G.; Yee, D.; Sorensen, L. B. Nature 1994, 368, 444-446. (45) Fenter, P.; Cheng, L.; Rihs, S.; Machesky, M.; Bedzyk, M. J.; Sturchio, N. C. J. Colloid Interface Sci. 2000, 225, 154-165.

NaCl Solution-NaCl Crystal Interface

Figure 1. Left panel: Sketch of the rectangular simulation box. Periodic boundary conditions (pbc) are applied to all three directions. The NaCl crystal (shaded) is sandwiched between solution phases (unshaded) and vice versa (due to the pbc). Right panel: Density Fbulk of aqueous NaCl solution at 298 K and 1 bar vs salt concentration cbulk: solid line, experimental values;56 plus signs, bulk solution at NaCl crystal with SPC/ E-P water model; diamonds, homogeneous solution with SPC/ E-P model; triangle, bulk solution at NaCl crystal with SPC/E model.

are aware of only two papers dealing with aqueous electrolytes at the (meta)stable NaCl(001) interface: purely classical MD simulation by Shinto et al.10 and a partially quantum-mechanical embedded cluster approach by Stefanovich et al.14 The authors do not use electrolytes at finite concentrations, but single sodium and chloride ions in aqueous phases at NaCl(001) interfaces. Water is treated either as a nonpolarizable sterical model10 or as a dielectric continuum.14 Therefore, we think simulations of moderately to highly concentrated aqueous NaCl solutions at the NaCl(001) interface, allowing the water molecules to respond to local electrical fields, can lead to new insights into structural and dynamical properties of electrolytes at uncharged ionic interfaces. The models simulated and methods employed in this study are specified in the following section 2. In section 3, the results, especially structural and dynamical properties of sodium and chloride ions in SPC/E-P water at the salt interface, are presented in detail. This section also contains a short comparison to a system with SPC/E water and dipole moment profiles of SPC/E-P water. Finally the results are discussed, and conclusions are summarized in section 4. 2. Models and Methods The heterogeneous systems considered here consist of crystalline sodium chloride in contact with aqueous sodium chloride solutions of varying concentrations or with neat water. The interfaces, lying perpendicular to the z-axis, are formed by NaCl(001) crystal layers in all cases. The simulation cell is depicted in the left panel of Figure 1. Since periodic boundary conditions are applied to all three directions, the crystal slab is surrounded by aqueous layers, and the aqueous layer is sandwiched between two crystal surfaces. Molecular Models. We mostly employ the polarizable extended simple point charge (SPC/E-P) water model by Sto¨ckelmann et al.,7 but the SPC/E model34 of Berendsen et al. (nonpolarizable) is also used. During simulations, the atomic charges of SPC/E-P water are adjusted to local electrical fields according to the fluctuating charge algorithm by Rick, Stuart, and Berne,46 which is based on the charge equilibration method by Rappe´ et al.47 Unlike the point charges of SPC/E water, the delocalized charges (46) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141-6156.

Langmuir, Vol. 18, No. 2, 2002 549

of the SPC/E-P model are characterized by atom-centered Slater-type orbitals. The nonpolarizable sodium and chloride ions possess a proton and electron charge, respectively. The distributions of ionic charges are described in terms of Slater-type orbitals such as those of the SPC/E-P water partial charges. The details of ionic interactions as well as Lennard-Jones, Born-Mayer, and electrostatic parameters for SPC/E-P water, ion-water, and ion-ion interactions are given elsewhere.7 But the reported Born-Mayer energy parameter for interactions between sodium and chloride ions has to be corrected, it should read: ANaCl ) 575400 kJ/mol. The same parameters are used for crystal and dissolved ions. The ionic interaction parameters fitted to SPC/E-P water by Sto¨ckelmann et al.7 are used throughout this work, even in the simulation employing SPC/E water. So the latter can serve only to elucidate how the water model affects the simulation results. Simulation Details. In our MD simulations, we use the AMBER 4.1 program package,48-50 treating long-range electrostatic interactions with particle mesh Ewald (PME) summation,51 extended by Sto¨ckelmann52 to include the fluctuating charge algorithm. In general, the temperature is kept at 298 K employing Berendsen’s thermostat,53 but a few runs at 323 and 348 K, respectively, are also reported (the time constant is τT ) 0.2 ps). The pressure (target value of 1 bar) is adjusted using Berendsen’s barostat53 (time constant τP ) 0.2 ps, isothermal compressibility κ ) 4.46 × 10-5 bar-1) by scaling the box lengths and particle positions anisotropically. The molecular geometry is fixed using the SHAKE algorithm.54 A time step of 1 fs is employed in all cases. A cutoff radius of 9 Å is applied to the direct sum part of PME summation and the LennardJones and Born-Mayer pair interactions. The crystal ions are not locally fixed and can thus vibrate and (in principle) dissolve from the crystal. If all dissolved ions are located near a single interface at the beginning of a run, equilibration can be monitored by comparing ion density profiles close to the lower and upper interface, respectively. Our simulations usually cover a period of 8-9 ns for time averaging after 3-4 ns of equilibration. An MD run of overall 12 ns takes about 40 days on an EV6 processor of a Compaq workstation XP1000 operating at a rate of 500 MHz. Simulated Systems. The MD runs carried out are given in Table 1. In general, a solution of NaCl in SPC/ E-P water in contact with a six-layered salt crystal is simulated at 298 K. But there is also a system (I) with neat water as liquid, a system (XIV) containing SPC/E water, a system (VIII) with a nine-layered crystal, and two runs (IXb, IXc) at elevated temperatures. The other systems differ from each other by the number of dissolved (47) Rappe´, A. K.; Goddard, W. A., III J. Phys. Chem. 1991, 95, 33583363. (48) Pearlman, D. A.; Case, D. A.; Caldwell, J. W.; Ross, W. S.; Cheatham, T. E., III; Ferguson, D. M.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman, P. A. AMBER 4.1; University of California, San Francisco, 1995. (49) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179-5197. (50) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1996, 118, 2309. (51) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089-10092. (52) Sto¨ckelmann, E. Ph.D. Thesis, Universita¨t Mainz, Germany, 1999. (53) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (54) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987.

550

Langmuir, Vol. 18, No. 2, 2002

Oyen and Hentschke

Table 1. List of MD Runsa MD run

water model

I II III IV V VI VII VIII IXa IXb IXc X XI XII XIII XIV

SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E-P SPC/E

no. of no. of no. of NaCl concn in temp crystal crystal dissolved bulk phase (K) ions layers ion pairs (mol/L) 298 298 298 298 298 298 298 298 298 323 348 298 298 298 298 298

864 600 384 600 384 864 600 900 384 384 384 600 600 600 600 600

6 6 6 6 6 6 6 9 6 6 6 6 6 6 6 6

0 5 5 10 10 20 20 20 20 20 20 40 80 160 320 40

0 (0.2-0.3) 0.25 0.55 0.55 (≈1) 1.02 1.11 1.07 1.09 1.05 2.16 4.17 7.72 12.8 2.34

aThe solution phase contains 1000 water molecules in all cases. The crystal layers perpendicular to the z-axis consist of square arrays of 64, 100, and 144 ions, respectively. Approximate bulk concentrations in parentheses indicate systems without bulk solution phases due to insufficient thicknesses of solution layers.

ion pairs and thus by salt concentration or by the number of ions per crystal layer, which takes on values of 64, 100, and 144, respectively. With growing number of ions per crystal layer, the cross-sectional area parallel to the interface is increased and therefore the statistics of ion density profiles etc. is improved, but the thickness of the solution layer is decreased because the number of water molecules is kept constant. In this case Lz (cf. Figure 1) may be too small to obtain a bulk phase region with constant and identical concentrations of cations and anions (runs II and VI). The NaCl(001) interfaces of all the systems listed in Table 1 turn out to be metastable during the simulation period: No ions are dissolved from the crystal into neat water or diluted solutions, and no crystallization occurs, not even for the most concentrated solution. The metastable crystals are barely affected by the solution phases: The density profiles of crystal sodium and chloride ions along the z-axis are virtually the same and only a little widened at the interface, regardless of the electrolyte concentration; the crystal layers are not permanently polarized. Since the experimental solubility of NaCl in water at 298 K lies at 5.4 mol/L,55 corresponding to 111 ion pairs in 1000 water molecules, the solutions of MD runs XII and XIII have to be regarded as oversaturated. So these systems are considered in special cases only. We note that the bulk solution densities Fbulk obtained at 298 K match the experimental densities at the corresponding bulk NaCl concentrations cbulk56 very well (cf. the right panel of Figure 1, where data of homogeneous systems consisting of 1000 SPC/E-P water molecules and 20 or 80 NaCl ion pairs, respectively, are also included). However, the simulated bulk phase densities Fbulk are underestimated at elevated temperatures (as in the case of neat SPC/E-P water52). Average Potentials of Mean Force. We define average potentials of mean force (PMFs) w(z) based on atomic density profiles g(z) in the following manner

w(z) ) -kBT ln g(z)

(1)

where kB is Boltzmann’s constant and T is the temperature. The density profiles are normalized to the bulk phase, i.e., g(z) ) 1 for large z. (All the average PMFs in this

paper are given as molar quantities, i.e., Boltzmann’s constant is substituted by the universal gas constant.) Degreve et al.20 calculate PMFs of ions in homogeneous solution from radial distribution functions in a way analogous to eq 1. The term PMF is chosen because the gradient of the PMF between two structureless (spherically symmetric) particles in homogeneous solution equals the mean force which the two particles exert on each other if they are locally fixed and the average is taken over configurations of all the other particles (cf. Chandler57). Considering such a dissolved structureless particle close to a corrugated interface, Chandler’s derivation can be transferred directly to heterogeneous systems only if the particle distribution function g(x,y,z) depends not only on the distance z from the interface but also on the position in the layer parallel to the interface (coordinates x, y). This leads to a relation analogous to eq 1 between the PMF w(x,y,z) and g(x,y,z). Since our average PMFs w(z) are obtained (according to eq 1) from density profiles g(z) averaged over the cross-sectional area A parallel to the interfaces, i.e., g(z) ) 1/A ∫A∫g(x,y,z) dx dy, they are related to the corresponding PMFs w(x,y,z) by the following equation:

w(z) ) -kBT ln

(

1 A

∫A∫ exp

)

-w(x,y,z) dx dy kBT

(2)

Some authors (e.g., Rose et al.58) call these average PMFs “free-energy profiles”. 3. Results 3.1. Structural and Dynamical Properties of Ions Dissolved in SPC/E-P Water. Ion Density Profiles and Average Potentials of Mean Force. Particle distributions perpendicular to the interface, i.e., along the z-axis, are characterized by atomic density profiles normalized to simulated bulk phase densities and nonnormalized concentration profiles, respectively. Here z ) 0 is defined as the mean z-position of the ions of the outermost crystal layer. (Note that the density profiles of crystal sodium and chloride ions are virtually identical.) Concentration profiles of sodium and chloride ions obtained from MD runs with bulk NaCl concentrations of ≈1 mol/L are displayed in Figure 2a. The sodium ion profiles are nearly identical, and the same applies to the chloride ion profiles, even in the case (run VI) in which there is no bulk phase with constant and identical sodium and chloride ion concentrations. This means, on the one hand, the concentration profiles are not affected significantly by the number of ions at the crystal surface, i.e., by the interfacial area and the thickness of the solution layer in the range considered here (runs VI, VII, IXa). Comparing systems VII and VIII, which differ by the number of crystal layers only (six and nine layers, respectively), on the other hand, shows that a six-layered crystal is sufficient to avoid finite-size effects. The most prominent features of the concentration profiles in Figure 2a are the intense selective adsorption of sodium ions and lack of chloride ions at the interface. For comparison, the normalized density profiles for oxygen and hydrogen atoms of water are displayed in Figure 2b. These profiles are (55) Stephen, H., Stephen, T., Eds. Solubilities of Inorganic and Organic Compounds; Macmillan: New York, 1963; Vol. 1, part 1, pp 108 + 109. (56) Lide, D. R., Ed. CRC Handbook of Chemistry and Physics, 81st ed.; CRC Press: Boca Raton, 2000; p 6-8. (57) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987; pp 201 + 202. (58) Rose, D. A.; Benjamin, I. J. Chem. Phys. 1991, 95, 6856-6865.

NaCl Solution-NaCl Crystal Interface

Figure 2. (a) Concentration profiles cNa(z) and cCl(z) of sodium and chloride ions along the z-axis perpendicular to the interface, obtained from MD runs with bulk NaCl concentrations of about 1 mol/L. (b) Densities gO and gH of oxygen and hydrogen atoms normalized to bulk phase densities (i.e., gObulk ) gHbulk ) 1) vs distance z from the interface (MD run IXa). (c) Normalized ionic density profiles gNa and gCl along the z-axis for MD run X. The error bars are based on six independent trajectories by splitting the 9 ns production run into three sections of 3 ns each and considering the lower and upper interfaces separately.

only little affected by dissolved ions, except in the cases of supersaturated bulk solutions (not shown here), and they exhibit the three hydration layers already observed by Sto¨ckelmann et al.7 The adsorbed sodium ions in the distance range from 1.5 to 3.5 Å with a maximum sodium concentration at 2.5 Å from the interface are located in the first hydration layer and between the first and second hydration layers, respectively; i.e., there can be no water molecules between the adsorbed and crystal ions. Starting configurations with no and with too many adsorbed ions at an interface, respectively, lead to the same ionic distribution; i.e., the reversibly adsorbed ions are in equilibrium with the bulk phase. Since the crystal layers do not exhibit any permanent electrical dipole moment in the z-direction, the sodium adsorption and lack of chloride at the interface cannot be caused by electrostatic interactions with a polarized crystal. The fine structure of the ionic concentration profiles is observed in the case of a ≈2 M bulk solution, i.e., in Figure 2c, more clearly than in Figure 2a: The normalized sodium ion density profile shows not only a pronounced peak at 2.6 Å but also an extremely low density at 3.5-4.5 Å, a small peak at 5.2-5.4 Å, and another minimum at 6.2 Å; there are two small, but significant, maxima in the chloride ion density profile at 4.6-5.0 Å and at 6.2-6.6 Å, i.e., in the range between 4 and 7 Å with higher than bulk phase chloride concentration. The three sections of the production run of 3 ns each, used for calculating the error bars in Figure 2c, are statistically independent, because the largest residence times found for low-mobility ions in the 4 Å thick solution layers adjacent to the interfaces are less than 1.5 ns (see below). Although there are 40 dissolved ion pairs in system X, only four to five sodium ions are adsorbed at each interface on the average over 9 ns, which explains why at least a ≈2 M bulk NaCl solution is required to verify the fine structure of ionic

Langmuir, Vol. 18, No. 2, 2002 551

Figure 3. (a) Density profiles gNa(z) of sodium ions normalized to bulk phase NaCl density along the z-axis. (b) Normalized densities gCl(z) of chloride ions vs distance z from interface. (c) Average potentials of mean force (PMFs) wNa(z) of sodium ions. (d) Average PMFs wCl(z) of chloride ions. The bulk NaCl concentrations have the following values: dotted line, 0.55 mol/L (MD run IV); solid line, 1.02 mol/L (run VII); dashed line, 2.16 mol/L (run X); long dashed line, 4.17 mol/L (run XI).

density profiles. This fine structure cannot be attributed to direct interactions between dissolved and crystal ions, so the hydration structure at the interface must strongly affect ionic distributions at least in the interfacial distance range between 4 and 7 Å. Normalized sodium and chloride ion density profiles, obtained for different bulk NaCl concentrations ranging from about 0.5 to 4 mol/L, are displayed in parts a and b of Figure 3. All the sodium ion density profiles exhibit an adsorption maximum at z ≈ 2.5 Å with 5- to 8-fold bulk density. These profiles match each other closely, even regarding the fine structure. The chloride ion density profile, on the other hand, depends on the bulk NaCl concentration (interplay of hydration effects and Coulombic forces as the probable mechanism for the different behavior of sodium and chloride ions discussed below). The broad first peak of the chloride ion density profile is widened with decreasing bulk concentration due to the increasing effective interaction range of the electrostatic forces. Moreover, the fine structure, clearly visible at a bulk concentration of ≈2 mol/L, is less obvious in cases of about 1 or 4 M bulk solutions, and there is no such fine structure at a bulk concentration of ≈0.5 mol/L. The low normalized density of chloride ions at a distance of ≈3 Å, i.e., directly at the interface, grows with increasing bulk concentration. Only in the case of ≈13 M (oversaturated) bulk solution, however, does this density reach values sufficiently large for determining ionic mobilities (see below). The density profiles of sodium ions at distances above ≈4 Å barely depend on bulk NaCl concentration; i.e., they are not significantly affected by interionic electrostatic forces but are determined by the hydration structure alone. Contrary to that, the chloride ion density profiles in this range are affected by Coulombic interactions also, leading to the observed concentration dependence (widening of the first peak with decreasing concentration). We cannot distinguish between hydration and electrostatic effects upon chloride distribution, but we

552

Langmuir, Vol. 18, No. 2, 2002

Oyen and Hentschke

mention that the characteristic lengths of both effects are of the same order of magnitude in our simulations: The hydration layers have an overall thickness of ≈5 Å, and the Debye length equals about 1.5-4 Å in 4 to 0.5 M (homogeneous) NaCl solution at 298 K. The Debye length rD of (1,1)-electrolytes is calculated by rD ) (RT/2F2cbulk)1/2, where  is the permittivity of water, R is the universal gas constant, T is the temperature, F is Faraday’s constant, and cbulk is the molar concentration.59 (Of course, these Debye lengths are merely crude estimates because DebyeHu¨ckel theory applies to electrolytes at concentrations below 0.001-0.01 mol/L only.) The average potentials of mean force for sodium and chloride ions, obtained from the density profiles in panels a and b of Figure 3 with eq 1, are displayed in panels c and d of Figure 3. The small absolute values of these PMFs are remarkable. For instance, the intense adsorption of sodium ions corresponds to a minimum PMF value of only -5 to -4 kJ/mol at an interfacial distance of ≈2.5 Å. However, this picture neglects the dependency of the potential of mean force on the x- and y-coordinates (see above). Moreover, the adsorption and desorption of ions can hardly be characterized without explicitly considering current (local) hydration structures. The integral over the difference between the sodium and chloride ion concentrations in the range from the interface to a specific interfacial distance z, i.e., ∫0z(cNa(ζ) - cCl(ζ)) dζ, which is proportional to the net ionic charge between the interface and a plane parallel to the interface at the distance z, is plotted versus the distance z in Figure 4. The range of positive values above ≈2 Å, caused by sodium ion adsorption and incomplete shielding by chloride ions, decreases with increasing concentration. While the curves for bulk NaCl concentrations up to 1 mol/L (runs II, IV, VII) exhibit continuous decline at distances greater than ≈3 Å, i.e., the distance of maximum net ionic charge, the curve for ≈4 M bulk solution (run XI) shows oscillatory behavior in this range, and the oscillations even grow with further increasing the bulk concentration above the experimental solubility of NaCl (not shown here). Depending on bulk concentration, the transitional range from the interfacial region to the bulk phase lies at interfacial distances between 10 and 20 Å in the cases considered here. Mobility of Ions. Ionic mobilities perpendicular to the interface are characterized by ionic residence times τ in layers parallel to the interfaces. Here τ is calculated via

N(t) ) N(0) exp

-t τ

(3)

where N(t) is the average number of ions remaining inside the layer (without excursion) after time t of originally N(0) ions. Note that the residence time τ is obtained from the linear range of the logarithmic plot of the data. Since, in general, residence times strongly depend on layer thickness, we employ 4 Å thick layers in all cases. Residence times of sodium and chloride ions, respectively, are plotted vs the mean distance of the layer from the closest interface in parts a and b of Figure 5. The residence times of adsorbed sodium ions range from 500 to 1100 ps. The other sodium ions exhibit residence times of only 3070 ps; i.e., the mobility of sodium ions is reduced by a factor of about 20 due to interfacial adsorption. Contrary to that, residence times of chloride ions vary by a factor of only 3 (20-60 ps). These results are barely affected by (59) Atkins, P. W. Physical Chemistry, 3rd ed.; Oxford University Press: Oxford, 1986.

Figure 4. Differences of sodium and chloride concentrations integrated in the range from 0 Å to z, i.e.,∫0z(cNa(ζ) - cCl(ζ)) dζ, vs distance z from interface.

Figure 5. (a) Average residence time τNa of sodium ions in 4 Å thick layers parallel to the interfaces vs mean distance z of the layers from the closest interface. Bulk concentrations cbulk range from about 1/4 to ≈1 mol/L. (b) Average residence time τCl of chloride ions in 4 Å thick layers parallel to the interfaces vs mean distance z of the layers from the closest interface. Systems, bulk concentrations, cbulk, and corresponding symbols are the same as those in (a). (c) Mean diffusion coefficient Dxy for lateral diffusion of sodium and chloride ions, respectively, in 4 Å thick layers parallel to the interfaces vs mean distance of the layers from the closest interface. Data for MD run VI only.

the bulk NaCl concentration in the range from about 1/4 to 1 mol/L (cf. Figure 5) and up to ≈4 mol/L (not shown). The diffusion coefficient Dxy of (two-dimensional) lateral ionic motion is determined via the Einstein relation54

Dxy )

〈[r b(t) - b r (0)]2〉 1 lim 4 tf∞ t

(4)

where t is the time and b r(0) and b r(t) are the particle positions projected onto the interfacial xy-plane at the time origin and at the time t, respectively, and averages are taken over time origins and all the sodium or chloride ions found in a 4 Å thick layer parallel to the interfaces. Dxy is obtained from the slope of the linear range in the plot of mean square projected distance 〈[r b(t) - b r(0)]2〉 vs time t. The resulting diffusion coefficients of sodium and chloride ions are plotted vs the mean layer distance from the closest interface in Figure 5c. Here we restrict our consideration to the system of MD run VI, but we mention

NaCl Solution-NaCl Crystal Interface

Langmuir, Vol. 18, No. 2, 2002 553

Figure 6. Average hydration number Nhydr of sodium ions in 1 Å thick layers parallel to the interfaces vs mean layer distance z from the closest interface at three different bulk NaCl concentrations cbulk.

that the lateral diffusion coefficients for interfacial distances below 8-10 Å do not depend on the extension in the z-direction of the solution layer. The diffusion coefficient of sodium ions is decreased by a factor of 5-10 due to adsorption at the interface, while the diffusion coefficients of chloride ions vary by a factor of only 2-3; i.e., adsorption reduces sodium ion mobility parallel to the interface also, but the effect is less pronounced than in case of mobility perpendicular to the interface. The density of chloride ions directly at the interface (at distances up to 3.5 Å) reaches values sufficient for detailed analysis in case of ≈13 M bulk solution (MD run XIII) only. Although bulk phase mobilities are lower in supersaturated solutions than in unsaturated ones, this allows us to compare adsorbed sodium ions with chloride ions in the same distance range (z e 4 Å). While chloride ions are more mobile perpendicular to the interface than sodium ions (residence times τ 200-300 ps for Cl-, ≈800 ps for Na+), the lateral diffusion coefficient Dxy of chloride ions (0.0004 Å2/ps) is much less than the one for sodium ions (0.01 Å2/ps). The mean distances of diffusion, d, calculated by the Einstein-Smoluchowski relation59 for diffusion in two dimensions, i.e., d ) (4Dxyτ)1/2, are therefore ≈0.6 Å for Cl- and ≈6 Å for Na+; i.e., adsorbed sodium ions diffuse parallel to the crystal surface, while chloride ions close to the interface merely vibrate. (For comparison, the mean distance between neighboring like-charged ions of the NaCl(001) surface equals 4.0 Å.) This different behavior of sodium and chloride ions can be explained by the hydration structure: The oxygen atoms of the water molecules in the first hydration layer are located close to the crystal sodium ions, thus forming a two-dimensional lattice.52 While the adsorbed sodium ions occupying interstitials of this lattice close to the crystal chloride ions (see below) can diffuse without disrupting the lattice, the lateral diffusion of chloride ions requires displacement of oxygen atoms at the crystal sodium ions. Hydration Numbers. We define the hydration number of a sodium ion as the number of oxygen atoms which are within the first hydration shell defined by the first minimum of the sodium-oxygen radial distribution function (RDF). Normalizing this RDF, we take into account that the oxygen atoms cannot penetrate the crystal phase. The resulting RDF, which is virtually independent of the distance between the dissolved ion and the closest interface, exhibits a maximum at a Na-O distance of 2.2 Å and a minimum at 2.9-3.1 Å. Since the RDF shows a pronounced first minimum, the first hydration shell of sodium ions is well defined. In Figure 6, the mean hydration numbers of sodium ions in 1 Å thick layers parallel to the interfaces are plotted vs the average layer distance from the closest interface. The hydration numbers of sodium ions adsorbed at the interface range from 4.2 to 4.4. In all other cases, this number is between 5.0 and

Figure 7. Lateral distribution of solution ions close to the interface: 2D RDF g2 vs interionic distance projected onto the xy-plane, rNa-Na (upper panel), and rNa-Cl (lower panel), respectively (system XI). The fat solid curves correspond to ions close to the interface, i.e., zNa e 4 Å and zCl e 3.5 Å. For comparison, the fat dashed curves show lateral distributions of bulk phase ions (averages over six layers, each of which is 2 Å thick), and the thin solid lines characterize the outermost crystal layers.

5.2. Thus the hydration shell of a sodium ion loses on the average about one water molecule upon adsorption. This supports the picture that there are no water molecules between the adsorbed and the crystal ions. The normalized chloride-oxygen RDF, which also is virtually independent of the interfacial distance of the dissolved ion, shows a maximum at a Cl-O distance of 3.1-3.2 Å and a minimum at 3.8-4.1 Å. This RDF exhibits a smooth first minimum; i.e., the first hydration shell of chloride ions is badly defined, and exact hydration numbers analogous to those of sodium ions are not obtained. Lateral Structure. The lateral distributions of solution ions in layers parallel to the interfaces are characterized by projecting the ion positions onto the xy-plane and calculating distance distribution functions of the projections. These distribution functions are calculated analogous to bulk RDFs. Thus we call them “two-dimensional radial distribution functions” (2D RDFs). Since the accuracy of 2D RDFs strongly depends on the number of ions, all the 2D RDFs presented in this paper are obtained from a system with ≈4 M bulk NaCl solution (MD run XI). Even in this case (with overall 80 dissolved ion pairs), the average number of adsorbed sodium ions (zNa e 4 Å) per interface is only about 8. The corresponding number of chloride ions (with zCl e 3.5 Å) is still smaller, i.e., about 0.3. In this case, a crystal layer consists of 100 ions, i.e., 50 ion pairs. The upper panel of Figure 7 shows not only the 2D RDF of adsorbed sodium ions but also the 2D RDFs of sodium ions in bulk phase layers and of sodium ions in the outermost crystal layer. The 2D RDF of adsorbed ions exhibits neither the broad peak around 3.5 Å of the bulk phase 2D RDF nor the narrow peak at 4 Å of the crystal 2D RDF. At distances above 5-6 Å, the 2D RDF of adsorbed sodium ions exhibits bulklike behavior, except for the small maxima at 8 and 12 Å (and at 6 Å), which correspond to maxima of the 2D RDF in the crystal. The 2D RDF between adsorbed sodium ions and chloride ions with zCl e 3.5 Å is displayed in the lower panel of Figure 7. Also shown are the 2D RDFs between sodium and

554

Langmuir, Vol. 18, No. 2, 2002

Oyen and Hentschke

Figure 8. Lateral distribution of solution ions at the interface relative to crystal surface ions: 2D RDF g2 vs interionic distance projected onto the xy-plane, rNa-Cl (upper panel), and rNa-Na and rCl-Cl (both lower panel), respectively (system XI). Key: fat solid lines, RDFs between adsorbed sodium ions (zNa e 4 Å) and crystal chloride/sodium ions; fat dashed lines, RDFs between chloride ions near interface (zCl e 3.5 Å) and crystal sodium/ chloride ions; thin solid lines, the corresponding functions calculated for adjacent crystal layers (for comparison).

chloride ions in the bulk phase and in the outermost crystal layer. Overall, the 2D RDF of solution ions at the interface is similar to the corresponding bulk 2D RDF. The 2D RDFs between solution ions directly at the interface and their respective counterions in the outermost crystal layer are displayed in the upper panel of Figure 8. For comparison, the 2D RDF of sodium and chloride ions calculated for two adjacent crystal layers is also shown. All these functions exhibit a high and narrow peak at zero projected distance; i.e., solution ions are preferentially adsorbed at crystal counterions. The dissolved ions are less localized than the crystal ions, and there is a clear difference between adsorbed sodium ions and chloride ions at the interface: The 2D RDF between adsorbed sodium ions and chloride ions of the outermost crystal layer indicates a significant number of ion pairs with projected distances around 2 Å (also pointing at lateral diffusion of adsorbed sodium ions), where the other two 2D RDFs show a wide (nearly complete) void (i.e., no lateral diffusion of solution chloride ions). Even with projected distances above 2-3 Å, the 2D RDF of ions from neighboring crystal layers is more similar to the 2D RDF between solution chloride ions and sodium ions of the outermost crystal layer than to the one between adsorbed sodium ions and chloride ions of the outermost crystal layer. The 2D RDFs between solution ions directly at the interface and their respective co-ions in the outermost crystal layer as well as the 2D RDF between like ions of two neighboring crystal layers are shown in the lower panel of Figure 8. Especially the pronounced peaks at 3 and 6.5 Å of the 2D RDFs of sodium (or chloride) ions in neighboring crystal layers are clearly visible in the 2D RDF between solution and crystal chloride ions but not in the one between adsorbed and crystal sodium ions. These results indicate that the few solution chloride ions with interfacial distances zCl e 3.5 Å better match the NaCl crystal structure than the many adsorbed sodium ions; i.e., the adsorption of sodium ions cannot simply be

Figure 9. Normalized density profiles g(z) vs distance z from the nearest interface at three different temperatures: (a) sodium; (b) chloride; (c) oxygen; (d) hydrogen. Panels e and f, average potentials of mean force of ions: wNa(z) (e) and wCl(z) (f).

considered as partial crystallization due to interactions between adsorbed and crystal ions alone. The results obtained from MD runs with oversaturated bulk NaCl solutions (not given here) are virtually identical. Temperature Dependency. Atomic density profiles and corresponding average potentials of mean force (PMFs), obtained from MD runs of the same system at three different temperatures ranging from 298 to 348 K, are shown in Figure 9. Selective interfacial adsorption of sodium ions is observed in all cases, but it decreases with increasing temperature and reaches about half its room temperature value at 348 K. In general, the normalized densities of sodium and chloride ions in Figure 9a,b are smoothed out with increasing temperature, and especially the first peaks decrease significantly. The first maxima in the density profiles of oxygen and hydrogen (panels c and d), on the other hand, are barely affected by changing the temperature. Therefore, the effect of temperature on interfacial adsorption of sodium ions cannot be attributed to a drastical change of the hydration structure, but the adsorption must be exothermic. The average PMFs of sodium and chloride ions, calculated from the density profiles in panels a and b of Figure 9 with eq 1 and displayed in panels e and f, exhibit a weak temperature dependence. For instance, the minimum PMF value of sodium ions at 2.5 Å is increased by 1.4 kJ/mol only (from -4.9 to -3.5 kJ/mol) due to a temperature rise of 50 K. 3.2. Comparison of Water Models. Results, obtained from two MD runs with different water models, the polarizable SPC/E-P model (run X) and nonpolarizable SPC/E water (run XIV), but with the same ion-ion and ion-water interactions as well as the same crystal size and numbers of dissolved ions and water molecules, are

NaCl Solution-NaCl Crystal Interface

Langmuir, Vol. 18, No. 2, 2002 555

Figure 11. Average dipole moment µ of SPC/E-P water in 0.2 Å thick layers parallel to the interfaces vs mean distance z of the layers from the closest interface. (a) Bulk concentration dependence of dipole moment profile: dotted line, neat water; solid line, ≈1 mol/L; dashed line, 2.16 mol/L; long dashed line, 4.17 mol/L. Data from MD runs I, VI, X, and XI (all runs at 298 K). (b) Temperature dependence from MD runs with about 1 M bulk solutions: solid line, 298 K; dotted line, 323 K; dashed line, 348 K (runs IXa, IXb, and IXc).

Figure 10. Normalized atomic density profiles g(z) (upper panel) and average potentials of mean force w(z) (lower panel) of sodium (gNa, wNa), chloride (gCl, wCl), oxygen (gO, wO) and hydrogen (gH, wH) vs interfacial distance z: solid lines, SPC/ E-P water (MD run X); fat dashed lines, SPC/E water (MD run XIV). PMFs are given in kJ/mol.

presented in Figure 10. The most important finding is that the strong selective adsorption of sodium ions, observed in the case of SPC/E-P water, does not occur in the case of SPC/E water. The density profile of sodium ions in SPC/E water exhibits fine structure including two moderately high peaks at interfacial distances of 4.0 and 5.2 Å and a shoulder at 2-3 Å. The differences between the chloride ion density profiles obtained with the two water models are less pronounced though significant. Even in the case of SPC/E water, the net ionic charge between the interface and xy-planes with interfacial distances z below 6-8 Å, calculated by integrating the difference of the sodium and chloride ion densities over the interfacial distance in the range from 0 Å to z, is positive but much smaller than in the case of the polarizable model; i.e., sodium ions are adsorbed preferentially in both cases, and the effect is strongly enhanced by explicitly considering water polarization. The density profiles of oxygen and hydrogen atoms of SPC/E-P and SPC/E water are mostly independent of the water model employed: Especially the first peak of the oxygen atom density profile at 2.0-2.2 Å and the first peak of the hydrogen atom density profile at 2.6-2.8 Å with a shoulder at ≈2 Å are nearly identical for both water models. However, the smooth second and third peaks of the oxygen density profiles exhibit observable differences (in the range from 3 to 6 Å). The average PMFs obtained from the density profiles in the upper panel of Figure 10 via eq 1 are displayed in the lower panel of the same figure. While the average PMFs for oxygen and hydrogen atoms, respectively, are virtually identical (except for the oxygen PMFs at 3-6 Å), the average PMFs for sodium and chloride ions, respectively, differ significantly but little from each other. Even the potential of sodium ions at an interfacial distance of 2.6 Å is decreased by ≈8 kJ/mol only when polarization is explicitly introduced into the water model.

The comparison of ionic distributions obtained with the two water models shows that the water molecules not only affect but even determine the distribution especially of sodium ions. The strong interfacial adsorption of sodium ions, occurring in SPC/E-P, but not in SPC/E water, cannot be attributed to direct interactions between dissolved and crystal ions but is caused by the structure of hydration layers so that little differences in the water structure, which are not clearly visible in the oxygen and hydrogen atom density profiles, can cause huge differences in ionic distributions. So the NaCl(001) interface affects ionic distributions via hydration layers formed at the interface. 3.3. Dipole Moment Profiles of SPC/E-P Water. Figure 11a shows dipole moment profiles of SPC/E-P water and three NaCl solutions of various bulk concentrations at the NaCl(001) interface. The average dipole moment of pure water equals about 2.70 D in the bulk phase but declines to 2.63 D directly at the interface, which is in accord with the findings of Sto¨ckelmann et al.7 Considering water molecules at interfacial distances z above ≈5 Å, i.e., molecules in the range from the third hydration layer to the bulk phase, the dipole moment decreases slightly with increasing NaCl concentration, and it reaches a value of 2.67 D in ≈4 M bulk solution. Contrary to that, the average dipole moment of water molecules at z e 2.4 Å, i.e., directly at the interface, rises with increasing bulk NaCl concentration (valid for unsaturated bulk solution). In the case of ≈4 M bulk solution, the average dipole moment at z e 2.0 Å (2.68 D) even slightly exceeds the bulk phase value (2.67 D). All the dipole moment profiles exhibit some fine structure: a distinct maximum of 2.702.71 D at distances of 3.2-3.6 Å in the second hydration layer (2.8-5.0 Å), at least a weak maximum at 5.6-6.0 Å in the third hydration layer (5.0 - ≈7 Å), and two minima (partly weak) roughly at the borderlines of the first and second layer and of the second and third layer, respectively. The decline in bulk phase dipole moment with increasing salt concentration can be explained in terms of ions perturbing the water structure. However, the average dipole moment of water in ≈4 M bulk solution is reduced by 0.03 D only as compared to neat water, while the dipole moment of SPC/E-P water vapor (1.86 D7) lies ≈0.8 D below the liquid-phase value. This interpretation is supported by the general decrease of dipole moment due

556

Langmuir, Vol. 18, No. 2, 2002

to rising temperature, which is exhibited by the dipole moment profiles displayed in Figure 11b. When the temperature is increased from 298 to 348 K, the gross structure of dipole moment profiles is preserved, but the fine structure becomes less clear. 4. Discussion and Conclusions The molecular dynamics simulations presented in this work show that all the employed models of aqueous NaCl solutions at the NaCl(001) interface exhibit at least a weak selective adsorption of sodium ions. In the case of the SPC/ E-P water model, this adsorption is strong and exothermic. This observation is in accord with the positive sign of the surface charge of NaCl particles in saturated salt solution, found experimentally by laser-Doppler electrophoresis.15 The conclusions of Ohtaki et al.8 from their simulation of a NaCl crystal dissolving in water (cf. section 1) indicate that the water molecules can cause the adsorption of sodium ions even if the sodium ions are hydrated stronger than the chloride ions. The densities of bulk NaCl solutions obtained from our simulations at 298 K and 1 bar are in very good agreement with experiment. Modeling the adsorption of ions and SPC/E water molecules at the clean NaCl(001) surface, Shinto et al.10 obtained an adsorption energy of about -20 kJ/mol for both sodium and chloride ion adsorption at their respective counterions. Their reported energy of water adsorption strongly depends on the adsorption site: ≈-10 kJ/mol at a surface chloride ion and -43 kJ/mol at a sodium ion.9 Sto¨ckelmann’s et al.7 SPC/E-P water is also preferentially adsorbed at sodium ions of such a NaCl(001) surface, where the authors determined an adsorption energy of -36 kJ/ mol. These numbers indicate that sodium ion adsorption is favored over water adsorption at chloride ions of the NaCl(001) surface and that water adsorption is favored over chloride ion adsorption at surface sodium ions. The numbers are also in accord with the hydration structure observed in simulations of bulk water at the NaCl(001) interface employing the SPC/E9 and SPC/E-P model,52 respectively. The oxygen atoms of the first hydration layer are located near the crystal sodium ions in both cases. In the case of bulk aqueous phases in contact with the NaCl(001) interface, however, the adsorption of ions is not analogous to the one at the clean surface. For instance, the potential of mean force (PMF) obtained by Shinto et al. for sodium and chloride ions in SPC/E water at the NaCl(001) interface is generally much greater than the potential energy of the ions at the corresponding distance

Oyen and Hentschke

from the clean surface.10 The PMF of adsorbed sodium ions is even positive, and the sodium ion PMF exhibits energy barriers of about 10 kJ/mol, not found in the potential energy curve of sodium ions at the clean surface, which clearly displays the effect of hydration layers. We note that the average PMFs of ions calculated in our work are generally lower than those of Shinto et al. (called “free energy profiles” by them). Shinto et al.10 obtained their PMFs by integrating the forces acting on solute ions located above their respective counterions, while we calculated our average PMFs from ion density profiles averaged over the cross sectional area parallel to the interfaces. Modeling the adsorption of sodium and chloride ions and of water molecules at the NaCl(001) interface in contact with a continuum water model and at the clean NaCl(001) surface, respectively, Stefanovich et al.14 obtained much smaller absolute values of adsorption energies in the bulk water case than in the other one. The authors concluded that water molecules adsorbed at the NaCl(001)-water interface are rather like hydrated molecules perturbed by the interface than like water molecules adsorbed at the clean surface and perturbed by the aqueous phase. Nevertheless, the differences between the SPC/E model and its polarizable version, the SPC/E-P model, offer the following interpretation for the preferential adsorption of sodium ions close to the crystal surface. Water bonds quite strongly to the sodium ions of the crystal trying to optimize its structure under this constraint. This is consistent with the low concentration of chloride ions in the range below 4 Å from the surface. The solution sodium ions trying to occupy the favorable chloride sites of the surface are likely to strain the water structure near the surface. The resulting stress is more easily relaxed by a water model which has the additional freedom to redistribute its charges. This argument also supports the idea that a polarizable model is the better model in a strongly heterogeneous environment. However, there are no quantitative experimental data for discriminating between the two water models. So further experiments and simulations of experimentally stable interfaces are needed for resolving this question. Acknowledgment. Financial support by the Deutsche Forschungsgemeinschaft is gratefully acknowledged. We thank Dr. E. Sto¨ckelmann for making available his modified version of the AMBER 4.1 program package. LA011269E