J. Phys. Chem. B 2006, 110, 20491-20501
20491
Molecular Dynamics Simulations of Electrolyte Solutions at the (100) Goethite Surface Sebastien Kerisit,*,† Eugene S. Ilton,† and Stephen C. Parker‡ Chemical Sciences DiVision, Pacific Northwest National Laboratory, Richland, Washington 99352, and Department of Chemistry, UniVersity of Bath, ClaVerton Down, Bath BA2 7AY, U.K. ReceiVed: June 12, 2006; In Final Form: August 5, 2006
Molecular dynamics simulations of electrolyte solutions in contact with a neutral (100) goethite (R-FeOOH) surface were used to probe the structure of the mineral-water interface and gain insight into the adsorption properties of monovalent ions. Three electrolyte solutions were considered: NaCl, CsCl, and CsF. The electrolyte ions were chosen to cover a range of ionic sizes and affinities for the aqueous phase. The molecular dynamics simulations indicate the presence of a structured interfacial region resulting from the strong interaction of water with the mineral surface. The specific arrangement and preferred orientation of water that arise from this interaction create adsorption sites in the interfacial region, i.e., as far as 15 Å away from the surface, and hence give rise to a strong correlation between the water and ion distributions. The structure of the hydrated ion, its effect on the water arrangement at the interface, and the strength of the ion-water bond are found to be key factors that determine the location and extent of ion adsorption at the interface. Additionally, in all simulations, we find a build up of positive charges near the surface due to cation adsorption, which is compensated by an accumulation of anions in the next few angstro¨ms. This creates an excess of negative charges, which is in turn compensated by an excess of positive charges, and so on. As we modeled a neutral surface, the structure of the electrolyte distribution arises from the complex interplay of the interactions between the surface, water, and the electrolyte ions rather than from the need to neutralize a surface charge. In addition, our simulations indicate that the electrolyte distribution does not resemble that of a classical electrical double layer. Indeed, our calculations predict the presence of several condensed layers and oscillations in the net charge away from the surface.
Introduction A molecular-level description of the mineral-water interface is essential to our understanding of many geochemical processes such as adsorption, the incorporation of trace metal impurities in crystals, and crystal growth and dissolution. A key feature of the mineral-water interface, and one that is incorporated into most phenomenological models of any process occurring at the mineral-water interface, is the distribution of electrolyte ions in the vicinity of the surface. The structure of the mineralsolution interface is usually based on the electrical double layer (EDL) model proposed by Gouy, Stern, and others (see ref 1 and references therein). In their model, the EDL is composed of the Stern (or condensed) layer where counterions are adsorbed to the mineral surface by electrostatic forces or by formation of complexes with surface atoms and the Gouy (or diffuse) layer where ions move freely and the net electric charge decreases exponentially with distance. Recent progress in experimental techniques and computer simulations has made it possible to evaluate the validity and limitations of these models through microscopic characterization of the EDL. For example, X-ray standing wave studies of the EDL at the rutile (110)-water interface2-6 have investigated the structure of the condensed layer in different conditions for a range of mono- and multivalent cations and have determined the partitioning of these cations between the condensed and diffuse layers. These techniques provide insight into the behavior of ions at mineral-water * Corresponding author. E-mail:
[email protected]. † Chemical Sciences Division, Pacific Northwest National Laboratory. ‡ Department of Chemistry, University of Bath.
interfaces but also have limitations. Direct measurements of low-Z elements are difficult, and density profiles in the vicinity of mineral surfaces have so far only been published for water.7,8 Therefore, molecular dynamics (MD) simulation is a valuable complementary approach that can provide a basis for interpretation of experimental data as well as a detailed molecular-based picture of the interface. For example, Rustad and co-workers9-12 have successfully modeled the protonation of surface sites at the mineral-water interface and Spohr13,14 has calculated the distribution of electrolyte ions in the vicinity of metal surfaces and the effect of surface charge on the structure of the EDL. Recently, we have used MD simulations based on an empirical potential model to study the interfacial region between both charged and uncharged iron (oxyhydr)oxide surfaces and a sodium chloride aqueous solution.15,16 Ultimately, our goal is to model the interaction of contaminants with the surfaces of such minerals in order to predict and assess the risks associated with heavy metal pollution of soils and sediments. Indeed, iron (oxyhydr)oxides have an excellent capacity to accumulate metal pollutants by adsorption17 and are therefore prime candidates for reducing heavy metal mobility. However, we first need to conduct a detailed investigation of the structure and dynamics of geochemically relevant mineral-water interfaces in the absence of contaminants. In this paper, we probe the behavior of electrolyte solutions at the goethite (R-FeOOH)-water interface as a function of cation and anion type. We consider electrolyte ions that span a wide range of sizes and affinities for the aqueous phase to evaluate the essential parameters that dictate the structure of the interface. Relatively high electrolyte concentrations are used to keep the length of the molecular
10.1021/jp0636569 CCC: $33.50 © 2006 American Chemical Society Published on Web 09/19/2006
20492 J. Phys. Chem. B, Vol. 110, No. 41, 2006 simulations reasonable and still obtain reliable statistics. Nonetheless, we go beyond our previous calculations by increasing the simulated time by 1 order of magnitude. The present work supports our previous calculations15,16 which suggested that the distribution of electrolyte ions is structured even when the surface charge is neutral and that the explicit treatment of water molecules is essential to correctly capture the effects of the mineral surface on the electrolyte solution. In addition, our simulations lead to new findings, namely, that the structure of the hydrated ion, solvent rearrangement upon adsorption, and the strength of the ion-water interaction are important factors in determining the adsorption behavior of the electrolyte ions and that the partitioning of these ions between the interfacial and bulk regions does not follow the trend expected from hydration enthalpies. Computational Methods MD simulations of a neutral (100) goethite slab in contact with NaCl, CsCl, and CsF electrolyte solutions were performed using the computer code DL_POLY.18 These calculations are based on the Born model of solids,19 in which the atoms of a system are represented as point-charge particles that interact via long-range Coulombic forces and short-range interactions. The latter are described by parametrized functions and represent the repulsion between electron-charge clouds, the van der Waals attraction forces, and, where appropriate, many-body terms such as bond bending. In addition, our model accounts for the polarizability of anions via a shell model first described by Dick and Overhauser,20 where the polarizable ions consist of two particles, a core and a shell, coupled by a harmonic spring. MD simulations of NaCl solutions in contact with a NaCl surface by Oyen and Hentschke21 have shown that their polarizable version of the SPC/E water model yields a different sodium distribution at the interface compared to the nonpolarizable SPC/E model. Although no experimental data are available for comparison, they suggested that the more sophisticated polarizable model might be required in heterogeneous systems. The potential model used in this work is the same as described previously15 except for the addition of potential parameters to model the interactions of fluoride and cesium ions with the remainder of the system. The core and shell charges and the core-shell spring constant of fluoride ions as well as the parameters for the short-range fluoride-fluoride interaction were taken from Catlow et al.22 The potential parameters describing the interaction between fluoride ions and the oxygen ions of the goethite mineral were derived by de Leeuw et al.23 based on the parameters for fluoride ions by Catlow et al.22 The potential models for the fluoride-water oxygen and fluoridewater hydrogen interactions were fitted to electronic structure calculations of small fluoride-water gas-phase clusters reported in the literature (see Table II of the Supporting Information). These potentials were then tested against experimental and theoretical data of the structure of the fluoride hydration shell in solution and were found to give good agreement (see Table III of the Supporting Information). Similarly, the potential parameters for the cesium-water interaction were derived to reproduce the structure and energetics of cesium-H2O gas-phase clusters and were subsequently tested against experimental and theoretical data on the structure and dynamics of the ion’s hydration shell in solution. Again, good agreement was obtained (see Tables IV, V, and VI of the Supporting Information). Both fluoride- and cesium-water models slightly overestimate the ion’s coordination number; however, this is in keeping with the relatively high calculated density of the water model used in
Kerisit et al.
Figure 1. Top view of the (100) goethite surface in a vacuum (left) from an energy minimization and in aqueous solution (right) from a snapshot of a MD simulation showing chemisorbed water molecules. Iron atoms are green, oxygen atoms red, and hydrogen atoms white.
this work.24 Indeed, due to the weak cesium-water interaction and hydrogen bonding between water molecules in the fluoride hydration shell, the water-water interaction plays an important role in determining the coordination number of these ions, unlike, for example, divalent cations where the ion-water interaction is much stronger and the orientation of the first shell water molecules prevents hydrogen bonding. Keeping the F parameter of the cesium-water oxygen Buckingham potential fixed (see Table I of the Supporting Information for detail of all the potential forms used in this work), we scaled the A parameter to the charges of the lattice and hydroxyl oxygen ions to obtain the cesium-lattice oxygen and cesium-hydroxyl oxygen potentials, following standard mixing rules. Similarly, the F parameter of the cesium-water oxygen Buckingham potential was kept fixed and the A parameter was fitted to the lattice constant of CsCl and CsF to obtain the cesium-chloride and cesium-fluoride potentials, respectively. These potentials were then evaluated against experimental physical constants of the two minerals and were found to give very good agreement (see Tables VII and VIII of the Supporting Information). Finally, the A parameter of the iron-fluoride potential was fitted to the hexagonal FeF3 crystal structure while using the same F parameter as the iron-lattice oxygen and hydroxyl oxygen potentials (see Table IX of the Supporting Information). All the potential parameters used in this work are given in Table I of the Supporting Information. The mineral slab, which contained 192 FeOOH units, had a thickness of approximately 20 Å and a surface area of 18.80 × 18.87 Å2 as previously.15 The (100) surface was exposed on both faces of the mineral slab (the referred space group is Pnma and, therefore, this surface corresponds to the (010) surface in space group Pbnm, which is also often used9,17). We only considered the hydroxyl-terminated surface as it is the most stable termination according to our static calculations.15 The (100) surface is one of the most stable low-index surfaces of goethite.25 In addition, as the most stable termination is fully hydroxylated, we did not have to assume a particular degree of hydroxylation, hence making the surface termination unambiguous. The terminal hydroxyl groups are coordinated to two cations from the Fe-O layer below, as shown in Figure 1. In a vacuum, the topmost iron atoms are in 5-fold coordination and thus the surface is corrugated with grooves between hydroxyl rows that are 4.7 Å wide. The mineral slab was put in contact with a water slab approximately 65 Å thick, which contained 960 water molecules. Therefore, the simulation cell dimensions were 18.80 Å × 18.87 Å × 84.00 Å and periodic boundary conditions were applied in three dimensions. Three different solution compositions were considered, namely, NaCl, CsCl, and CsF solutions. In each case, 16 cations and 16 anions were introduced randomly in the water slab, which corresponds to a salt concentration of
Electrolyte Solutions on Goethite
J. Phys. Chem. B, Vol. 110, No. 41, 2006 20493
1.2 mol‚dm-3, i.e., about twice that in seawater. The electrolyte concentration is much lower than the solubility of any of the three compounds (6.2, 11.3, and 37.7 mol‚dm-3 at 25 °C for NaCl, CsCl, and CsF, respectively26). The atomic velocities were initially scaled for 5 ps to the target temperature. Next, the system was allowed to equilibrate for 200 ps before the production run started and data were collected for 10 ns, which corresponds to 50 × 106 MD steps. We monitored the projection on the normal to the surface of the trajectories of randomly selected sodium, chloride, and water oxygen ions and found that these ions explored all regions of the aqueous solution within the time scale of the simulation (see Figure I of the Supporting Information). This suggests that the length of the simulations reported here is sufficient for the systems to reach equilibrium and to minimize kinetics effects on the ions’ distribution. Furthermore, as will be discussed in detail later, the dynamics of water and ion exchange between the surface and the solution are at least 1 to 2 orders of magnitude shorter than the simulated time. All the simulations were performed in the NVT ensemble (constant number of particles, constant volume, and constant temperature) at 300 K and zero applied pressure. The temperature was kept constant by used of the Nose´-Hoover thermostat.27 The electrostatic interactions were calculated by means of the Ewald summation method28 with a real space cutoff of 8 Å. The same cutoff was used for the shortrange interactions. Interaction between periodic images in the direction normal to the surface is a possible issue; however, Spohr29 showed that in simulation cells of sufficient size the use of three-dimensional rather than two-dimensional periodic boundary conditions does not significantly affect calculated properties of the interface and aqueous solution, except for integrated properties such as the electrostatic potential, which can show quantitative discrepancies. Additionally, in previous work,16 we applied the correction of Yeh and Berkowitz30 to similar simulation cells and found the correction to be small and it was thus deemed unnecessary. The particles’ trajectories were generated by the Verlet leapfrog algorithm with a time step of 0.2 fs. The shells were given a small mass of 0.2 au, and their motion was treated as that of the core following the adiabatic shell model first introduced by Mitchell and Fincham.31 All profiles normal to the surface such as density profiles were obtained using a bin height of 0.1 Å and radial distribution functions were calculated with a bin size of 0.01 Å. Results and Discussion Electrolyte Distribution at the Interface. MD simulations of the NaCl, CsCl, and CsF electrolyte solutions in contact with the (100) goethite surface were run for 10 ns. Figure 2 shows the ion distributions and water densities normal to the surface obtained from averaging the data from both faces of the mineral slab and normalizing to the concentration of a homogeneous solution. The water densities were obtained from the position of water oxygen atoms. The surface (distance ) 0 Å) is defined as the plane passing through the center of the uppermost layer of iron atoms. Figure 2 indicates the presence of two regions in the solution. The first region, which will be referred to as the layered or structured region thereafter, extends to about 15 Å and shows oscillations in water and ion densities. It is then followed by a 20 Å thick region where the density of each species is nearly constant and equal to that in a homogeneous solution. The damped oscillatory behavior of the water density is in very good agreement with previous simulations of this mineral15 and others32,33 as well as experimental results.7,8 The oscillations in the ions’ distribution correlate strongly with the
Figure 2. Ion distributions and water density as a function of distance from the surface in the simulations of NaCl (top), CsCl (middle), and CsF (bottom). Average of both faces is shown.
fluctuations in the water density. This can be explained by the three-dimensional structure and specific orientation of water at the interface which create preferred locations in the layered region for the ions to reside and hence give rise to layered ion distributions. Therefore, analogy with the EDL model suggests that the condensed layer can span several hydration layers and hence consists of several adsorption sites. It is worth mentioning that cesium displays the least amount of correlation beyond the first adsorption peak, in accord with its structure breaking properties. Due to the weak interaction between cesium ions and water, little water reorientation around cesium ions is required as they enter the layered region and hence rapid exchange between the different adsorption sites is expected. Figure 2 also suggests that the presence of electrolyte ions has very little effect on the water density fluctuations in the direction normal to the surface.
20494 J. Phys. Chem. B, Vol. 110, No. 41, 2006
Kerisit et al.
TABLE 1: Percentage of Electrolyte Ions Adsorbed in the First Density Peak and within 10 Å from the Surface in the Three Simulations species
simulation
first density peak
within 10 Å
Na+ ClCs+ ClCs+ F-
NaCl NaCl CsCl CsCl CsF CsF
16.5 7.9 8.5 4.9 6.5 1.4
27.9 28.4 20.1 24.3 16.8 17.2
The most striking result from Figure 2 is that sodium displays a much stronger adsorption peak at the interface than cesium. The proportion of ions associated with the first density peak and within the layered region was calculated from the integration of the distribution curves in Figure 2 and is reported in Table 1. On the basis of our calculated hydration enthalpy, -366 and -156 kJ/mol for sodium and cesium, respectively (see Table X of the Supporting Information for details of these calculations), one would expect sodium to have a greater affinity for the aqueous phase than cesium; however, since the surface is terminated by hydroxyl ions and chemisorbed water molecules (Figure 1), sodium is more strongly bound at the interface than cesium. Moreover, due to its small size, sodium can adsorb in a site near the surface with only minimum solvent rearrangement (see next section) and hence sodium adsorption is energetically favored relative to cesium. Also, due to its smaller size and shorter bonds with water and hydroxyl oxygen atoms, sodium adsorbs closer to the surface than cesium (first adsorption peak centered at 3.40 Å for sodium in the NaCl simulation and at 4.70 and 4.90 Å for cesium in the CsCl and CsF simulations, respectively). We find that, generally, cations adsorb closer to the surface than anions, except in the simulation of CsF; however, in this simulation, cesium is in excess in the third hydration layer and thus all the simulations show a cation excess in the first few angstro¨ms from the surface. The nature of the anion slightly affects the position and the size of the first cesium peak (position, 4.70 and 4.90 Å for the CsCl and CsF simulations, respectively; size, 2.12 and 1.49 for the CsCl and CsF simulations, respectively). In contrast to the cations, fluoride ions remain in the bulk solution to a larger extent than chloride ions in accord with their calculated hydration enthalpy: -676 and -533 kJ/mol, respectively. The difference can be explained in part by the fact that the adsorption of chloride ions in the layered region forces less disruption of the water structure than the adsorption of the smaller fluoride ions (see next section). Comparing the anion distributions obtained in the CsCl and CsF simulations, we find that fluoride ions, due to their smaller size, adsorb closer to the surface than chloride ions (first adsorption peak centered at 5.10 and 4.40 Å for chloride and fluoride ions, respectively). In addition, we find that the chloride density in the first 10 Å is greater than that of fluoride ions (24 and 17%, respectively). The anion distributions obtained in the CsCl and NaCl simulations indicate that the nature of the cation clearly influences the anion distributions. Indeed, to compensate for the larger excess of positive charges at the interface in the NaCl simulation, a greater density of chloride ions is present in the layered region with the difference in density being mostly located in the first adsorption peak. For each simulation, the charge due to the ions is integrated to obtain the electric field created by the electrolyte solution in the direction normal to the surface according to Poisson’s equation (Figure 3, top). The charge density due to the ions is
Figure 3. Electric field resulting from the distribution of the electrolyte ions for the three solutions (top) and that resulting from water in the NaCl and CsCl simulations together with the electric field due to NaCl for comparison (bottom).
defined as
Fq(z) )
∑i qiFi(z)
(1)
where Fi(z) is the density of ion i at height z and qi its charge. The electric field can be obtained from the charge density
E(z) )
1 -0
∫zz Fq(z′) dz′
(2)
0
where 0 is the permittivity of vacuum. The three curves show the same feature, namely, that the excess of cation density near the surface is compensated by an excess of anions in the next layer; however, this region extends beyond the first point of neutrality, which occurs between 7 and 10 Å, and thus results in a slight overall negative charge. This is in turn compensated by an excess of positive charges. This phenomenon gives rise to oscillations in the integrated charge. The period and the magnitude of these oscillations differ to a great extent between the three curves. Indeed, there is a 2 orders of magnitude difference between the electric field generated at the interface with the NaCl electrolyte solution and that with the CsCl and CsF electrolyte solutions. The behavior of the electric field at the interface indicates the presence of a response of the electrolyte ions to the structure of the mineral-water interface; however, it should be noted that, as the surface modeled in this work is neutral, the presence of such a response does not originate from the need to neutralize the surface charge but rather from the complex interplay of the mineral-water, mineral-electrolyte, and electrolyte-water interactions. In
Electrolyte Solutions on Goethite
J. Phys. Chem. B, Vol. 110, No. 41, 2006 20495
addition, the shape of these curves differs from the expected exponential decay of the EDL’s diffuse phase. This could be due to the high electrolyte concentrations employed in this work. Indeed, MD simulations of NaCl electrolyte solutions in contact with the NaCl (001) surface showed a transition from exponential decay to oscillations of the electric field at high salt concentrations.21 Similarly, MD simulations of molten and gaseous sodium chloride by Keblinski et al.34 showed that “strongly coupled regimes” develop such oscillations. Water is expected to make a nonuniform contribution to the overall electric field due to the specific orientation of water molecules at the interface and their dipolar nature. Therefore, we calculated the electric field considering all the water molecules using eqs 1 and 2, and we found that water generates high-frequency oscillations in the electric field over the first 15 Å from the surface (Figure 3, bottom). In addition, the amplitude of the oscillations is much greater than that generated by the electrolytes alone. Overall, the total electric field at the interface is dominated by the contribution from water, with slight modifications attributable to NaCl and CsCl. In conclusion, our results highlight the importance of the contribution of water to the properties of the interface and strongly suggest that an explicit representation of solvent molecules is required in microscopic models of the mineralsolution interface. In addition, our simulations indicate that the distribution of electrolyte ions cannot be directly rationalized from their affinity for the aqueous phase. Having described the partitioning of the electrolyte ions in the different regions of the interface, we now focus on the adsorption of these ions at the mineral surface. Adsorption at the Interface. Before addressing the adsorption of ions at the surface, we need to consider the structure and dynamics of the goethite-water interface. As mentioned previously, surface iron atoms are in 5-fold coordination and are thus undercoordinated in a vacuum (Figure 1). When the surface is in contact with an aqueous solution, water molecules adsorb between surface hydroxyl rows and complete the octahedral coordination shell of surface iron atoms (Figure 1). First, to investigate the structure of the layered region, we calculated the two-dimensional surface densities of the uppermost iron atoms, the surface hydroxyl ions, and the first five hydration layers (Figure 4). The iron-water bond is strong, which localized the chemisorbed water molecules and inhibited exchange with other water molecules during the simulations. However, because our water model is not spontaneously dissociable, the potential dissociation of chemisorbed water molecules is beyond the scope of this paper. Water molecules are adsorbed in well-localized sites even beyond the first hydration layer. In the fourth and fifth layers (i.e., between about 6 and 9 Å from the surface), the water density is converging to a uniform distribution in the direction of the surface hydroxyl rows but some localization remains in the direction perpendicular to these rows. Second, to investigate the dynamics of water at the interface, we computed the residence time of water in each hydration layer, as shown in Table 3. Residence times were obtained from the residence-time correlation function of water in each layer, which is defined as32,35
〈R(t)〉 )
〈
1
Nt
∑ θi(0)θt(t)
N0 i)1
〉
(3)
where N is the number of water molecules in a particular hydration layer and θi(t) is the Heaviside function, which is 1 if the ith water molecule is in the hydration layer at time t and
Figure 4. Two-dimensional densities from the simulation of NaCl: in all three plots, surface iron (green); top, surface hydroxyl oxygen (blue) and chemisorbed water oxygen (gray); middle, second (orange) and third (gray) hydration layers; bottom, fourth (orange) and fifth (gray) hydration layers. Color intensity increases with density.
0 otherwise. A water molecule was counted as having left the hydration layer if it had done so for a least 2 ps.35 This is to allow for water molecules which temporarily leave the hydration layer and come back without entering the remainder of the solution for a significant amount of time to be counted in the residence time. The residence time can then be obtained by integration of 〈R(t)〉
τ)
∫0∞ 〈R(t)〉 dt
(4)
20496 J. Phys. Chem. B, Vol. 110, No. 41, 2006
Kerisit et al.
TABLE 2: Comparison of the Ion-Oxygen Radial Distribution Functions Obtained When the Ions Are in Bulk Water and When Adsorbed (from the simulation of CsCl and NaCl for the cations, from the simulation of CsF and NaCl for the anions, and from the simulation of CsCl for oxygen)a sodium ligand hydroxyls 1st layer 2nd layer 3rd layer 4th layer 5th layer total bulk
cesium
1st max
1st min
N
2.73 2.30 2.33 2.31
3.13 3.50 3.20 3.17
0.84 1.96 1.88 2.10
2nd 1st max max
1st min
3.48 5.03 4.03 5.03
3.53 3.93 4.47 4.09 4.56 4.35
3.18 3.18 3.16 3.16 3.20 3.19
n
0.14 1.36 3.95 1.77 3.36 1.37 5.94 11.81 2.29 3.19 6.23 4.33 3.18 4.08 10.52
fluoride 2nd 1st max max
1st min
n
4.58 6.23 5.38 5.38 5.18 6.13
3.24 3.63 3.41 3.31
1.09 3.54 1.42 2.42
2.63 2.65 2.65 2.66
chloride 2nd 1st max max 4.23 4.26 4.78 4.50 4.51
3.16 3.15 3.15 3.15 3.15
1st min 3.85 4.41 3.88 4.29 4.27
n
0.91 4.00 1.90 3.90 1.86 8.47 12.57 5.88 2.63 3.28 8.10 4.63 3.13 3.93 12.10
3rd layer water oxygen 2nd 1st max max
1st min
4.92 4.62 5.41 5.49 5.36 6.11
3.92 4.32 3.97 4.38 4.30
2.97 2.94 2.97 2.96 3.00
n
1.27 4.05 1.70 3.77 1.37 12.16 5.38 2.93 4.18 12.77
2nd max 4.73 6.10 5.36 5.50 5.43 5.87 5.68
a The surface RDF is divided between oxygen atoms of the different hydration layers and those of surface hydroxyls. The table shows the first and second maxima and first minimum in the ion-oxygen radial distribution functions together with coordination number at the first minimum (see Figure 2 in the Supporting Information for the full RDF plots).
TABLE 3: Water, Anion, and Cation Residence Times (ps) in Each Hydration Layer (average of both faces is shown)a species
simulation
1
2
3
4
5
6
H 2O H 2O H 2O FClClCs+ Cs+ Na+
CsF CsCl NaCl CsF CsCl NaCl CsF CsCl NaCl
40000 65000 22000
177 317 384
23 51 69 99 111 202 57 98
16 24 23 48 60 74
6 9 9 20 22 29
9 10 8 20 12 12
423
a
Numbers in italics are not converged and are predicted from exponential fitting of the residence time correlation function.
Due to the limited number of desorption events of chemisorbed water molecules, the residence-time correlation function of water in the first hydration layer is not fully converged and the residence time is thus obtained by fitting the correlation function to an exponential function.35 There is some indication that sodium adsorption at the interface enhances the exchange rate of chemisorbed water molecules although longer simulations and/or alternative methods would be required to confirm this result. Comparing the residence time of water in the second and third hydration layers with that of the electrolyte ions (Table 3) indicates that water molecules are interacting quite strongly with the surface at these distances and hence compete with ions for adsorption sites. The residence times of water in the second and third hydration layers are longer when chloride ions are present rather than fluoride ions and when sodium ions are present rather than cesium ions. This correlates with the fact that chloride and sodium concentrations are greater at the interface than that of fluoride and cesium ions, respectively. There is no significant difference between the three simulations beyond the third hydration layer. Finally, to investigate further the dynamics of the interface, we computed the diffusion coefficient of water and the electrolyte ions across the water slab from the mean-square displacement (MSD) of these species (Figure 5)
1 D(z)t ) 〈|ri(t) - ri(0)|2〉 6
(5)
where D(z) is the self-diffusion coefficient of the species of interest at height z above the surface, t is the time, and ri(t) is the position of the species at time t. As water and the electrolyte ions explore regions with different diffusion coefficients within the time scale of the simulations, we only considered short intervals, i.e., 5 ps, as introduced by Marrink and Berendsen.36
Figure 5. Diffusion coefficient of water and electrolyte ions as a function of distance from the surface in the simulation of NaCl, together with water density.
The geometric center of the 5 ps trajectory determined the position along the normal to the surface. Figure 5 shows a clear correlation between water density and diffusion coefficient, i.e., water diffusion is high in region of low water density and vice versa. The electrolyte ions display the same correlation between the density and diffusion profiles. Oscillations in the diffusion coefficient of water and the electrolyte ions disappear from about 15 Å away from the surface. We suggest that this distance can be seen as the boundary between the bulk solution and the interface for the surface and conditions modeled in our simulations. Indeed, beyond this point, the structure and dynamics of the aqueous solution can be deemed identical to those of the bulk. Water diffusion coefficient converges to 2.6 × 10-9 m2/s in good agreement with the diffusion coefficient obtained from a bulk water simulation32 (2.3 × 10-9 m2/s) and from experiment37 (2.3 × 10-9 m2/s). Sodium and chloride diffusion coefficients converge to 1.8 and 1.7 × 10-9 m2/s, respectively, in accord with bulk diffusion coefficient (1.4 and 1.3 × 10-9 m2/s, respectively) The CsCl and CsF simulations gave almost identical results for water and diffusion coefficients of 2.2 and 1.1 × 10-9 m2/s for cesium and fluoride ions, respectively (data not shown). To investigate the adsorption behavior of the two cations, we computed two-dimensional surface densities similar to those shown for water where only the cations associated with the first adsorption peak were considered (Figure 6). Analysis of Figure 6 reveals that cesium ions exhibit two adsorption sites: site A, above a surface chemisorbed water molecule, and site B, above
Electrolyte Solutions on Goethite
Figure 6. Two-dimensional densities from the simulations of NaCl (top), CsCl (middle), and CsF (bottom): surface iron (green), surface hydroxyl oxygen (blue), chemisorbed water oxygen (gray), cesium in first density peak (orange), and sodium in first density peak (orange). Cesium adsorption sites A and B are also shown (middle). Color intensity increases with density.
an iron row, midway between two iron atoms and close to a surface hydroxyl ion. The former is more populated. In contrast, sodium ions only show one adsorption site: above a chemisorbed water row, halfway between two chemisorbed water molecules. Another important difference between the two cations is that sodium adsorption sites appear well-separated, whereas cesium adsorption sites are not, indicating exchanges between these sites. This suggests that lateral diffusion of cesium ions
J. Phys. Chem. B, Vol. 110, No. 41, 2006 20497 on the surface could take place whereas sodium ions have to first desorb from the surface to diffuse laterally. To evaluate the strength of the cation-surface interaction, we calculated the residence time of the two ions when associated with the first density peak (Table 3). First, we find that the nature of the anion can affect the residence time (cesium residence time is 57 and 98 ps in the CsF and CsCl simulations, respectively). Additionally, sodium is calculated to remain adsorbed for as long as half a nanosecond or about 1 order of magnitude longer than cesium, hence demonstrating the formation of a much stronger interaction with the surface. To explain the differences reported so far between the two cations, we propose that two of the main factors controlling the adsorption behavior of the electrolyte ions are the extent of water rearrangement due to the ion’s hydration shell structure and the strength of the ionsurface interaction relative to the water-surface interaction. As mentioned above, because the surface is terminated by hydroxyl ions and water molecules, we expect ions with a high hydration enthalpy to interact strongly with the surface. However, adsorption of these ions might not be energetically favored if it requires significant disruption of the water hydrogen bond network and adsorption of ions that reduce the net number of hydrogen bonds will have a smaller driving force. Indeed, using vertical positions from Figure 2 and lateral positions from Figure 4 and Figure 6, it can be shown that sodium does not displace any water molecule upon adsorption but instead, due to its small hydrated ionic radius, adsorbs at an “interstitial” site in the water structure where little or no water reorientation is required and as a result does not significantly reduce the number of hydrogen bonds. In contrast, cesium occupies a water oxygen site when adsorbed in site A, which requires displacing a water molecule and reducing the number of hydrogen bonds. Cesium cannot occupy a sodium interstitial site as its hydrated ionic radius is too large. This suggests that the strength of the cesium-surface interaction is too weak to compensate for the large solvent rearrangement that would be required for cesium to adsorb at this site. To evaluate whether the sodium site is not suitable for cesium or whether it was not accessed in our simulations due to a large barrier to adsorption, we replaced adsorbed sodium ions by cesium ions in the NaCl simulation and resumed the MD simulation. We found that all cesium ions rapidly diffused away from the sodium site ruling out any limiting kinetics effect on the cesium distributions. A related issue is the surface effect on the cations’ hydration shell. The cation-oxygen radial distribution functions (RDFs) obtained for cations in bulk water and those associated with the first density peak are shown in Table 2 (see also Figure II of the Supporting Information for the RDF curves). We find that sodium is on average coordinated by two water molecules in each of the first, second, and third hydration layers. Its coordination is only slightly reduced from that in the bulk (from 6.23 to 5.94). In addition, the RDF indicates a weak interaction with surface hydroxyl groups as shown by a broad peak centered at 2.73 Å (Figure II of the Supporting Information). Although sodium seems to retain most of its hydration shell upon adsorption, the chemisorbed water molecules only rarely exchange with the bulk solution on the time scale of the simulation, and thus sodium must replace part of its hydration shell by chemisorbed water molecules to adsorb at the surface. Cesium does not appear to be as strongly bound to the chemisorbed water molecules as sodium does. However, there is a well-resolved peak in the cesium-surface hydroxyl RDF at the bulk Cs-O distance of 3.18 Å (Figure II of the Supporting Information). We can conclude, from Figure 6, that this peak
20498 J. Phys. Chem. B, Vol. 110, No. 41, 2006 arises from bonding with the surface when cesium ions are adsorbed at site B. The fact that this peak integrates to 0.14 hydroxyl oxygens gives us an estimation of the partitioning between the two adsorption sites. To adsorb in site B cesium ions have to displace a water molecule in the second hydration layer. Table 3 shows that water molecules in this layer have long residence times, which limit the density of cesium at this site. For both cations, the position of the first maximum in the RDF is not significantly affected when the cations are adsorbed at the interface whereas the position of the second peak varies greatly. This confirms results obtained previously32 on the adsorption of metal ions at the calcite-water interface, namely, that there is competition between the surface and the adsorbed cation for water molecules. Generally, the cation-water interaction in the first shell is strong enough not to be affected by the surface whereas interactions with water molecules that belong to the second shell are much weaker and the second shell is distorted due to the structure imposed by the surface on water. Nonetheless, in the case of cesium, although the most probable distance remains largely unaffected, the first hydration shell is more distorted than that of sodium as can be seen from the position of the first RDF minimum. We now turn to the adsorption of the two anions. Again, we start by computing two-dimensional densities, only considering the anions associated with the first adsorption peak (Figure 7). As already mentioned, the adsorption of anions at this height above the surface is less extensive than that of the cations except for chloride when sodium is present. Also, we have already discussed how, when sodium is present rather than cesium, the chloride density at the surface is increased; however, Figure 7 indicates that the increased density is not due to a difference in the number of adsorption sites. Rather, the residence time of chloride ions doubles in the simulation of NaCl compared to that of CsCl (Table 3). In addition, Figure 7 shows that chloride and fluoride ions adsorb at the same site and that this site is also cesium’s adsorption site A, referred to as site A hereafter. Therefore, both of the anions, cesium, and water compete for site A. Visualization of the MD trajectories revealed that at any point of the simulation most chemisorbed water molecules form hydrogen bonds with surface hydroxyls and thus their oxygen ions are oriented toward the solution, as shown in Figure 1. Therefore, little reorientation of these molecules is required upon cesium adsorption whereas, when both fluoride and chloride ions occupy site A, chemisorbed water molecules rotate to form a new hydrogen bond with the adsorbed anions. Table 2 also shows data from the water oxygen RDFs with the different hydration layers when a water molecule occupy site A. The fact that the oxygen-first layer water RDF is similar to that of cesium but not to that of fluoride or chloride confirms that there is minimal reorientation of the chemisorbed water molecules upon cesium adsorption. As for the cations, the first maximum in the anion-water oxygen RDFs is only slightly affected; however, the positions of the first minimum and second maximum vary greatly. Again, this is strong evidence for distortion of the first and second hydration shells due to the presence of the surface. Table 2 allows us to characterize the structure of site A when cesium, fluoride, chloride, or water is occupying the site. Cesium, chloride, and water show very similar first and second shell structures indicating that, when cesium or chloride ions occupy site A in place of water, the disruption of the water structure is minimal. In contrast, the structure of the fluoride hydration shell induces large distortions
Kerisit et al.
Figure 7. Two-dimensional densities in the simulations of NaCl (top), CsCl (middle), and CsF (bottom): surface iron (green), surface hydroxyl oxygen (blue), chemisorbed water oxygen (gray), chloride in first density peak (orange), and fluoride in first density peak (orange). Color intensity increases with density.
in the water structure, which, we suggest, makes an endothermic contribution to the free energy of adsorption and limits fluoride adsorption. Finally, having discussed the influence of the surface on the ion-water interaction, we investigate how the cation-anion interaction differs at the interface from that in the bulk. Figure 8 shows the RDFs of the three cation-anion pairs for cations associated with the first density peak and in the bulk. In the aqueous phase, all three pairs form solvent-shared ion pairs
Electrolyte Solutions on Goethite
J. Phys. Chem. B, Vol. 110, No. 41, 2006 20499
Figure 8. Cation-anion radial distribution functions, whereby, the RDF was determined considering either all cations in the simulation or only those adsorbed at the surface (RDF labeled “surface”).
(SSIPs), in which water molecules bridge the two ions. However, only NaCl forms a significant amount of contact ion pairs (CIPs), in which there is no intervening water molecule between the two ions. The sodium-chloride RDF for bulk solution shows a main peak at 5.2 Å with a shoulder peak at 4.6 Å indicating that there are two possible configurations for the SSIP complex. Visualization of the MD trajectories indicates that the latter occurs when two water molecules are part of both hydration shells but do not lie directly between the two ions and the former corresponds to the situation where a single water molecule is positioned directly between the two ions, which results in a slightly longer Na+-Cl- distance. Both SSIP complexes show a greater density than the CIP complex, which indicates that they are lower in free energy than the CIP complex (see eq 6 below for the relationship between density and free energy). This is in very good agreement with many published MD simulations of NaCl pair formation.38-41 Interestingly, the three ion pairs show very different responses to the interfacial region. The amount of CsCl SSIP complexes is not affected by the surface, whereas that of CsF decreased by about half. NaCl shows an increase in the density of CIP complexes at the interface relative to the bulk. Furthermore, the SSIP peak increases at the interface whereas the shoulder peak decreases. This suggests that the interface structure does not allow for sodium-chloride ion pairs with two bridging water molecules to form but facilitates the formation of SSIP complexes where only one water molecule is shared. In summary, our calculations indicate that two of the key factors controlling the ion adsorption are the strength of the water-ion bond and the effect of the adsorbed ion on the water structure due to its size and the structure of its hydration shell. Finally, we observed that cation-anion pairs are affected differently by the surface depending on their structure relative to the topology of the surface and the water structure at the interface. In the next section, we investigate other possible ion pairs and how they are perturbed by the surface. Cesium-Cesium and Fluoride-Fluoride Pair Formation. Cation-cation and anion-anion RDFs show well-defined peaks in the cesium-cesium and fluoride-fluoride RDFs at 5.05 and 4.15 Å, respectively (Figure 9). This indicates the presence of free energy minima at these distances, which are consistent with configurations where two ions share water molecules of their first hydration shell. Therefore, our calculations suggest that ion pairs, and possibility larger clusters, are likely to form in aqueous solution. Figure 10 shows examples of these ion pairs from snapshots of the MD simulations. In the cesium pair,
Figure 9. Cation-cation (top) and anion-anion (bottom) radial distribution functions.
Figure 10. Snapshots of the cesium and fluoride ion pairs with their first hydration shell from the CsCl and CsF simulations, respectively.
several water molecules bridge the two cations through oxygen atoms whereas in the fluoride pair the bridging water molecules form hydrogen bonds with both anions. Such anion-anion and cation-cation pairs have some precedent in the literature. For example, MD simulations of fluoride-fluoride association in water that used an entirely different set of potential parameters42 showed a free energy minimum in the fluoride-fluoride potential of mean force (PMF) at a distance of 4.3 Å in good agreement with our results. It is important to note that the PMF can be obtained from the RDF using the relationship
W(r) ) -kBT ln(g(r))
(6)
where kB is Boltzmann’s constant, T is the temperature, and g(r) is the RDF density. It can be seen from eq 6 that a maximum in the PMF is a minimum in the RDF and vice versa. The same study42 reported an energy barrier to dissociation of 1.5 kcal/mol, which is consistent with our calculated value of 1.8 kcal/mol (see Figure III of the Supporting Information). To the best of our knowledge, cesium-cesium pairing in aqueous solution has not been investigated in the literature; however, Sutton and Sposito43 observed, in their MD simulations, sharing of water molecules by cesium ions adsorbed at the surface of Cs-smectite. Our calculations estimate the dissociation barrier
20500 J. Phys. Chem. B, Vol. 110, No. 41, 2006
Kerisit et al. TABLE 4: Cluster Size Distribution in Percentage of Total Number of Clusters Formed during Each Simulation (from the simulations of CsCl and CsF) cluster size
cesium (%)
fluoride (%)
2 3 4 5 6 7
70.86 22.12 5.41 1.41 0.16 0.03
91.22 8.67 0.12 0.00 0.00 0.00
consistent with their similar free energies of association (see Figure III of the Supporting Information). Finally, both cluster densities sharply decrease nearing the surface. However, fluoride shows a resurgence of cluster density in the adsorbed layer. The size distributions of fluoride and cesium clusters are shown in Table 4 and reveal that fluoride ions predominantly form pairs as well as a small number of three-membered clusters. The large majority of cesium clusters are pairs but with a greater proportion of three-membered clusters than fluoride. Cesium ions also form a small percentage of large clusters. To summarize, our MD simulations indicate the formation of cesium and fluoride pairs, which have different structural and dynamical properties. Fluoride ions form relatively tightly bound pairs bridged by hydrogen bonds whereas cesium ions can form larger clusters that are short-lived. Both ion pair cluster densities are influenced by the layered region. Conclusions
Figure 11. Excess cluster density as defined in the text as a function of distance away from the surface together with the ion and water distributions for reference.
of cesium pairs to be 0.8 kcal/mol. In addition, there is almost no energy barrier for association whereas fluoride ions show a barrier to association of about 1.2 kcal/mol. Therefore, cesium pairs and clusters are expected to form more readily but also to be much more short-lived than the fluoride pairs. To investigate the clustering process further, we calculated the cluster size distribution and the percentage of ions involved in a cluster as a function of distance from the surface. Ions were considered to be paired together if their interatomic distance was less than the first minimum in their respective RDFs. Clusters were identified by searching for unbroken chains of ion pairs. In Figure 11 the excess density of ion clusters is plotted as a function of distance from the surface. The excess density at a height z is calculated as
S(z) )
〈nClust. ions (z)〉 〈nTotal ions (z)〉 nhomo.
(7)
. where 〈nClust ions (z)〉 is the number of ions at height z above the surface that belong to a cluster, 〈nTotal ions (z)〉 is the total number of ions at height z, and nhomo. is the integration of the first RDF peak assuming a perfectly homogeneous solution. The term nhomo. normalizes the cluster density to the short ion-ion distances that simply occur from random fluctuations. Figure 11 shows several interesting features. First, the cesium curve is smoother than that of fluoride which is consistent with the fact that cesium clusters form more readily and also have shorter lifetime than fluoride clusters (see previous discussion). Second, in bulk solution, both ions form clusters in an excess of 2,
We performed MD simulations of three electrolyte solutions (NaCl, CsCl, CsF) in contact with a neutral (100) goethite surface. Our calculations indicate the presence of a structured interfacial region in the first 15 Å. In this interfacial region, the arrangement and orientation of water molecules are dictated by the structure and nature of the mineral surface, which also affects the diffusion properties of water to a great extent. This creates preferred locations for the ions to diffuse to and reside and thus the calculated ion distributions are directly correlated to the fluctuations in water density. In addition, electrolyte ions entering the interfacial region have to compete with the surface for water molecules and often with water molecules for adsorption sites. Our simulations indicate that, at the interface, the electric field due to the electrolyte ions shows an oscillatory behavior. In contrast, the presence of a structured electrolyte distribution at a neutral surface is not predicted by the classical EDL model. In addition, previous work15 showed that a small change in surface charge did not affect the shape of the electrolyte distribution, only the magnitude of the oscillations. Therefore, our calculations suggest that electrolyte solutions at neutral and charged surfaces show characteristics that differ in parts from the classical EDL model. Sodium and cesium are found to exhibit significantly different adsorption behaviors. Due to its smaller size, sodium occupies an interstitial site in the interfacial water structure that lies close to the surface. In contrast, cesium has to occupy a water adsorption site. This leads to a greater driving force for sodium adsorption as it does not disrupt the hydrogen bond network of interfacial water molecules. Similarly, fluoride ions, which greatly disrupt the interfacial water structure upon adsorption, show relatively small adsorption in the layered region in comparison to chloride ions which require less water rearrangement. Finally, we found the presence of a variety of ion pairs in the bulk of the aqueous solution. Cation-anion pairs tend to
Electrolyte Solutions on Goethite form solvent-shared complexes and no significant amount of contact-ion pair complexes except for sodium-chloride pairs. In addition, cation-cation and anion-anion ion pairs are observed for cesium and fluoride ions, respectively. Cesium pairs were found to form readily but to be short-lived, whereas fluoride pairs showed a higher barrier to association and to dissociation. The ion pairs were found to exhibit different behaviors in the interfacial region. We suggest that ion pairs that cause large distortions of the water structure will be thermodynamically unfavored at the interface. In the future, we wish to extend this work to the study of the adsorption of contaminants at the mineral-water interface and the effect of the interfacial structure on their adsorption properties. In addition, we plan to investigate a range of different minerals to evaluate to what extent the concepts put forward in this work depend on the nature of the mineral surface. Acknowledgment. The authors wish to acknowledge support for this work from the Pacific Northwest National Laboratory. This research was performed in part using the Molecular Science Computing Facility (MSCF) in the William R. Wiley Environmental Molecular Sciences Laboratory, a national scientific user facility sponsored by the U.S. Department of Energy’s Office of Biological and Environmental Research and located at the Pacific Northwest National Laboratory. Pacific Northwest is operated for the Department of Energy by Battelle. Supporting Information Available: A table summarizing all the potential parameters used in this work, eight tables comparing the structural parameters of small fluoride-water gas-phase clusters, the structure and dynamics of the first hydration shell of fluoride ions in aqueous solution, the energy and the cesium-water oxygen distance of cesium-water gasphase clusters, the structure and dynamics of the first hydration shell of cesium ions in aqueous solution, the physical constants of the minerals CsF and CsCl, and the lattice parameters of the mineral FeF3 obtained in this work with data published in the literature, a figure showing the projection on the normal to the surface of one randomly selected sodium, chloride, and water oxygen ion in the simulation of NaCl electrolyte solution, a table comparing the experimental and calculated enthalpies of hydration of the four electrolyte ions, a figure of the ion-oxygen radial distribution functions obtained when the ions are in bulk water and when adsorbed at the interface, and a figure showing the fluoride-fluoride and cesium-cesium potential of mean force. The material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Parsons, R. Chem. ReV. 1990, 90, 813. (2) Fenter, P.; Cheng, L.; Rihs, S.; Machesky, M. L.; Bedzyk, M. J.; Sturchio, N. C. J. Colloid Interface Sci. 2000, 225, 154. (3) Predota, M.; Bandura, A. V.; Cummings, P. T.; Kubicki, J. D.; Wesolowski, D. J.; Chialvo, A. A.; Machesky, M. L. J. Phys. Chem. B 2004, 108, 12049. (4) Predota, M.; Zhang, Z.; Fenter, P.; Wesolowski, D. J.; Cummings, P. T. J. Phys. Chem. B 2004, 108, 12061.
J. Phys. Chem. B, Vol. 110, No. 41, 2006 20501 (5) Zhang, Z.; Fenter, P.; Cheng, L.; Sturchio, N. C.; Bedzyk, M. J.; Predota, M.; Bandura, A.; Kubicki, J. D.; Lvov, S. N.; Cummings, P. T.; Chialvo, A. A.; Ridley, M. K.; Benezeth, P.; Anovitz, L.; Palmer, D. A.; Machesky, M. L.; Wesolowski, D. J. Langmuir 2004, 20, 4954. (6) Zhang, Z.; Fenter, P.; Cheng, L.; Sturchio, N. C.; Bedzyk, M. J.; Machesky, M. L.; Anovitz, L. M.; Wesolowski, D. J. J. Colloid Interface Sci. 2006, 295, 50. (7) Cheng, L.; Fenter, P.; Nagy, K. L.; Schlegel, M. L.; Sturchio, N. C. Phys. ReV. Lett. 2001, 87, 156103. (8) Fenter, P.; Cheng, L.; Park, C.; Zhang, H.; Sturchio, N. C. Geochim. Cosmochim. Acta 2003, 67, 4267. (9) Rustad, J. R.; Felmy, A. R.; Hay, B. P. Geochim. Cosmochim. Acta 1996, 60, 1563. (10) Felmy, A. R.; Rustad, J. R. Geochim. Cosmochim. Acta 1998, 62, 25. (11) Wasserman, E.; Rustad, J. R.; Felmy, A. R. Surf. Sci. 1999, 424, 19. (12) Rustad, J. R.; Wasserman, E.; Felmy, A. R. Surf. Sci. 1999, 424, 28. (13) Spohr, E. Electrochim. Acta 1999, 44, 1697. (14) Spohr, E. Solid State Ionics 2002, 150, 1. (15) Kerisit, S.; Cooke, D. J.; Marmier, A.; Parker, S. C. Chem. Commun. 2005, 24, 3027. (16) Spagnoli, D.; Cooke, D. J.; Kerisit, S.; Parker, S. C. J. Mater. Chem. 2006, 16, 1997. (17) Randall, S. R.; Sherman, D. M.; Ragnarsdottir, K. V.; Collins, C. R. Geochim. Cosmochim. Acta 1999, 63, 2971. (18) Smith, W.; Forester, T. R. DL_POLY is a package of molecular simulation routines, copyright The Council for the Central Laboratory of the Research Councils, Daresbury Laboratory at Daresbury, Nr. Warrington, 1996. (19) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; Oxford University Press: Oxford, U.K., 1954. (20) Dick, B. G.; Overhauser, A. W. Phys. ReV. 1958, 112, 90. (21) Oyen, E.; Hentschke, R. Langmuir 2002, 18, 547. (22) Catlow, C. R. A.; Diller, K. M.; Norgett, M. J. J. Phys. C.: Solid State Phys. 1977, 10, 1395. (23) de Leeuw, N. H.; Parker, S. C.; Hanumantha Rao, K. Langmuir 1998, 14, 5900. (24) de Leeuw, N. H.; Parker, S. C. Phys. ReV. B 1998, 58, 13901. (25) Cornell, R. M.; Posner, A. M.; Quirk, J. P. J. Inorg. Nucl. Chem. 1974, 36, 1937. (26) Physical Constants of Inorganic Compounds. In CRC Handbook of Chemistry and Physics, Internet Version 2006; Lide, D. R., Ed.; Taylor and Francis: Boca Raton, FL, 2006. (27) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (28) Ewald, P. P. Ann. Phys. 1921, 64, 253. (29) Spohr, E. J. Chem. Phys. 1997, 107, 6342. (30) Yeh, I. C.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155. (31) Mitchell, P. J.; Fincham, D. J. Phys.: Condens. Matter 1993, 5, 1031. (32) Kerisit, S.; Parker, S. C. J. Am. Chem. Soc. 2004, 126, 10152. (33) Kerisit, S.; Cooke, D. J.; Spagnoli, D.; Parker, S. C. J. Mater. Chem. 2005, 15, 1454. (34) Keblinski, P.; Eggebrecht, J.; Wolf, D.; Phillpot, S. R. J. Chem. Phys. 2000, 113, 282. (35) Impey, R. W.; Madden, P. A.; McDonald, I. R. J. Phys. Chem. B 1983, 87, 5071. (36) Marrink, S. J.; Berendsen, J. C. J. Phys. Chem. 1994, 98, 4155. (37) Sorenson, J. M.; Hura, G.; Glaeser, R. M.; Head-Gordon, T. J. Chem. Phys. 2000, 113, 9149. (38) Karim, O. A.; McCammon, J. A. J. Am. Chem. Soc. 1986, 108, 1762. (39) Rey, R.; Guardia, E. J. Phys. Chem. 1992, 96, 4712. (40) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757. (41) Geissler, P. L.; Dellago, C.; Chandler, D. J. Phys. Chem. B 1999, 103, 3706. (42) Dang, L. X. Chem. Phys. Lett. 1992, 200, 21. (43) Sutton, R.; Sposito, G. Geochem. Trans. 2002, 3, 73.
