Ion Distribution and Hydration Structure at Solid-Liquid Interface between NaCl Crystal and Its Solution.

The interface structure between NaCl crystal and its solution has been investigated at the saturated concentration of 298 K by molecular dynamics simulations. We have found that there are many fine structures at this complex interface. Near the surface of crystal, most of Na+ only coordinate with water molecules, while almost all Cl- coordinate with Na+ in addition to water molecules. An ion coordinating with more water molecules is farther away from the epitaxial position of lattice. As approaching to the interface, the first hydration shell of ions has the tendency of being ordered, while the orientation of dipole of water molecules in the first hydration shell becomes more disordered than that in the solution. Generally, the first hydration shell of Na+ is less affected by nearest Cl-, whereas the first hydration shell of Cl- is significantly affected by nearest Na+.

Introduction

In nature and daily life, the solid–liquid interface (SLI) between ionic crystals and solutions is ubiquitous. Many important physical phenomena and processes occur around SLIs, which are the key to understand the properties of two-phase boundary,[1−4] flotation,[5−8] crystal dissolution,[9−12] growth,[9,12−15] etc. For a long time, scientists have made great efforts to investigate the intrinsic feature of SLIs by various experimental techniques and computer simulations.

Due to the strong influence of the periodic character of crystal and long-range Coulomb interaction, the SLI between the ionic crystal and solution is usually more complicated than that of simple metal solid–liquid interfaces.[1,3,4,7,16,17] Compared to experimental methods, computer simulation is a powerful tool to obtain atomic scale information. The previous simulations have made important progresses for understanding the atomic scale structure of interfaces. It was found that both ions and water molecules exhibit a layerlike structure in the direction perpendicular to the interface.[15,17−31] Besides, the water molecules[19−21,24−26,28,30,32] and ions[14,19,20,24,25,28,30,32,33] at the interface are laterally ordered with respect to the crystal structure. The dipole moment of water molecules at the SLI also tends to be ordered.[17−19,21,24−26,34] Interestingly, ions at an SLI have much complexed adsorption structures,[18,20−24,27−29,31−33,35−37] which depends on the ions’ radius,[19,23,25,27,29,31−33] valence,[29,31] as well as the crystal surface structure.[21,27,29,30,38]

Hydration structures are the core structural unit for understanding the structure and physical properties of solutions. On the one hand, it is believed that the hydration shell still exists near the interface between crystal and solutions.[11,12,38−41] On the other hand, as approaching to interfaces, water molecules may rearrange with certain order.[7,17,34,42] It is not difficult to imagine that the hydration structure at SLI will be more complex under the two competing factors, namely, the tendency of maintaining hydration shell and the ordering tendency of water molecules. In addition to the impact on the structure, the ionic hydration shell is also the key to understand the transport of ions.[9,11,12,43] However, the ionic hydration structure itself has not been investigated systematically.

Interface between the NaCl and its solution is a typical example of this kind of SLI, which serves as a common system for studying complex interfaces.[4,7,11−14,17,34,42,44,45] In this work, the interfacial structure between NaCl crystal and saturated NaCl solution has been studied by molecular dynamics method. Our study has focused on the hydration structure of ions, especially on how the crystal surface affects both ions distribution and hydration structure and how the distribution of ions and ionic hydration structure interplay with each other. We find that there are many fine structures at the interface. The distribution of water molecules in the ionic hydration structure around the interface is strongly modified by the local environment of ions.

Results and Discussion

The density distribution of ions and oxygens along the direction is shown in Figure . Large density oscillations and layerlike structure are observed near the interface, which is universal in SLI systems. According to the density distribution of Na+ and Cl–, at least three layers can be identified. For the convenience of the following discussion, we assume that the first layer refers to the one closest to the crystal surface and next is the second and third layer in turn. In the first layer, the number of Na+ is obviously more than that of Cl–, which is consistent with previous work.[17] As demonstrated by Yang et. al.,[13] Na+ is easier to be adsorbed on crystal surface by replacing water molecules than Cl–.

Figure 1

Density distribution of H (green dot-dash line), O (red dash line), Na+ (blue dot line), and Cl– (magenta solid line) along the direction perpendicular to the interface, where the vertical dash line marks the epitaxial lattice position.

In the first layer, the density distribution of Na+ and Cl– with different nw’s is presented in Figure a,b, respectively. One can see that, in the first layer, Na+ prefers to coordinate with water molecules, while Cl– prefers to coordinate with Na+. Our calculations show that, around 39.7% of Na+ coordinate with Cl– and water molecules at the same time, while the rest Na+ do not coordinate with Cl– but only with water molecules (nc = 0). However, around 92.2% of Cl– coordinate with Na+ and water molecules at the same time and 7.8% of Cl– coordinate only with water molecules (nc = 0). Evidently, Na+ has stronger tendency to keep the hydration structure, which is consistent with previous results.[46]

Figure 2

Distribution of density of Na+ (a) and Cl (b) near the interface with various nw and nc, where the vertical dash line marks the epitaxial lattice position. nw denotes the number of water molecules in the first hydration shell, and nc is the number of the nearest counter ions.

The first peak of the density distribution for both Na+ and Cl– does not locate at the epitaxial lattice position (marked as the vertical dash line in Figures and 2). We find that the shift of the first peak is closely related to the local environment of ions, i.e., the coordination number of water molecules or counter ions around an ion. From Figure a,b, one can see that for both Na+ and Cl–, the less water molecules in the first hydration shell (FHS) are, the closer to the epitaxial lattice position is. If a Cl– is coordinated with more water molecules, it will be further away from the crystal surface. By contrast, a Na+ will be closer to the crystal surface if it is coordinated with more water molecules. The results can be roughly understood as follows: the distribution of ions in the first layer is mainly determined by two factors: (1) maintaining the ionic hydration shell and (2) driving ions to the ideal lattice position. Because of the mismatch between hydrate shell of ions and lattice structure, the former makes ions deviate from the ideal lattice, while the latter drives ions back to the ideal lattice. Therefore, with the decrease of nw, the ions will approach to the ideal lattice position gradually. In principle, the structure of interface is the result of competition between the energy and entropy. Our results show that the energy seems to be superior to entropy at the interface, indicated by the increased ordering tendency of both ions and water molecules.

From Figure , one can find another interesting feature, namely, the first layer is composed of a few sublayers. Each sublayer can be distinguished by the different nc or nw. We find that the sum of nc and nw is approximately a constant, which can be seen from Table . For Na+ and Cl–, the sum is around 5 and 6, respectively. In other words, if an ion coordinates with more counter ions, it will coordinate with less water molecules. If the density profile of the first layer is approximated by the Gaussian distribution, we find that the major difference among sublayers is the position of peak, not the half-width.

Table 1

Number of Water Molecules in the First Hydration Shell (nw) for Ions with Different Coordination Number of Counter Ions (nc)

	nc = 0	nc = 1	nc = 2	nc = 3
Na+	5.0 ± 0.05	4.0 ± 0.07	3.0 ± 0.13	1.9 ± 0.23
Cl–	6.0 ± 0.48	5.2 ± 0.32	4.3 ± 0.29	3.3 ± 0.69

The ion-oxygen radial distribution functions (RDFs) in the first three layers are shown in Figure a,b. For comparison, the ion-oxygen RDF in the solution is also shown. The position of the first peak in RDFs corresponds to the average radius of the FHS, and the area included by the first peak is proportional to nw. For both cation and anion, the interface has little effect on the average radius of the FHS, which can be seen in the inset of Figure a,b. For Na+ and Cl– in the first layer, as shown by the green line in Figure a,b, the first peak of ion-oxygen RDF is obviously lower than others, which indicates that the interface has an unneglectable effect on nw. However, this effect only works in the first layer, and for the second and third layers, the first peak of ion-oxygen RDF is approximately equal to that in solutions.

Figure 3

Ion-O RDFs. (a) Na+-O RDF in the first three layers and solution. (b) Same as (a) but for Cl–. (c) Na+-O RDF for two kinds of Na+ in the first layer, namely, coordinating with (labeled as type I) and without (labeled as type II) counter ions. Here, “total” denotes the sum of type I and type II. (d) Same as (c) but for Cl– in the first layer. The insets show the focus view of the first peak of ion-O RDF.

As mentioned above, ions with different nc values may have different local structures. To further clarify this difference, we divide both Na+ and Cl– in the first layer into two classes, namely, coordinating with (hereafter labeled as type I, namely, nc > 0) and without (hereafter labeled as type II, namely, nc = 0) counter ions. As shown in Figure c, the Na+-O RDFs for both type I and type II are similar, and the main difference appears in the first peaks. However, as shown in Figure d, there are some remarkable difference in the Cl–-O RDF between type I and type II. First, the first peak of Cl–-O RDF of type II is much higher than that of type I, due to more water molecules around type II Cl–. Second, near the first minimum of the Cl–-O RDF of type II, a small peak appears in the Cl–-O RDF of type I. These results may be the explanation why the radius of FHS of Cl– is different before and after dissolution.[47] Third, around the second peak of the Cl–-O RDF of type II, a minimum and a small peak appear in the RDF of type I.

The angle distribution of water molecules is the most intuitive physical quantity to describe the orientation of water molecules in the FHS. The distribution of cos(α) and cos(γ) near the interface is depicted in Figure . As expected, the water molecules uniformly distribute around ions in the solution and have no preferential orientation. However, as approaching to the interface, obviously the preferential orientation of water molecules in FHS appears. In other words, the FHS of ions tends to be ordered as approaching to the interface. When a Na+ approaches to the crystal surface, there are obviously two preferred directions, around cos(γ) = 0 and 1, and α prefers to stay around cos(α)= 0 and ±1. For Cl–, the distribution of cos(γ) is similar to that for Na+. On the contrary, the distribution of cos(α) for Cl– seems to have five preferred values of 0.0, ±0.71, and ±1.0, as approaching to the crystal surface.

Figure 4

Angle distribution for water molecules in FHS. (a, c) Distribution of cos(γ) and cos(α) of Na+. (b, d) Same as (a) and (c), respectively, but for Cl–.

The discussion above implies that the water molecules have certain ordering in the FHS. This kind of order is also closely related to the local coordination environment of ions. Figure presents the distribution of cos(α) and cos(γ) in the FHS for Na+ (a and b) and Cl– (c and d) in the first layer. The distribution of cos(α) and cos(γ) is similar for type I and type II Na+. The peak of type II Na+ at cos(α)=0 is sharper than that of type I, and the distribution of cos(γ) for type II just slightly shifts to the right relative to that of type I. However, coordinating with/without Na+ has a significant effect on the orientation of water molecules in the FHS of Cl–. The distribution of both cos(α) and cos(γ) for type I Cl– has more and sharper peaks than those for type II Cl–, which indicates that the water around type I Cl– is more ordered.

Figure 5

Distribution of cos(α) (a) and cos(γ) (b) in FHS of two types of Na+ in the first layer. (c, d) Same as (a) and (b), respectively, but for Cl–.

Dipole orientation of water molecules is another important physical quantity to describe the ordering of water molecules in the FHS. The dipole orientation of water molecules is characterized by φ angles (see Figure ), and the distribution of cos(φ) in the first layer is shown in Figure . The dipole of water molecules around both Na+ and Cl– in the solution prefers to point to a special direction, which is consistent with the previous results.[48,49] This special orientation, which minimizes the Coulomb interaction between the central ion and nearest water molecules, can be maintained in the interface beyond the first layer. However, in the first layer, due to the strong effect of crystal surface, the distribution of the dipole is no longer confined to a sharp region and is broadened. Namely, the dipole orientation of water molecules in the FHS is more disordered than that in the solution. For both type I and type II Na+, the distribution of cos(φ) in the first layer is similar, but different from that in solutions: a new peak appears near cos(φ) = 0.78, as shown in Figure a. For type II Cl–, the distribution of cos(φ) is similar to that in solution except a slight shift for the peak, as shown in Figure b. However, for type I Cl–, the distribution of cos(φ) is broadened remarkably relative to that of type II Cl–, and even some φ angles locate at the region where the interaction between Cl– and dipole becomes unstable locally.

Figure 8

Schematic plot of FHS of an ion, in which the purple, red, and blue spheres represent the central ions, oxygen, and hydrogen atoms, respectively; 1 is a unit vector from the central ion to an oxygen atom; 2 is another unit vector along the water molecule dipole; the , , and axes are defined in Figure . Rc is the cutoff radius of FHS.

Figure 6

Distributions of cos(φ) in the FHS of the first layer (black dot and red dash line) and of solutions (green solid line) for Na+ (a) and Cl– (b).

The change in the four angles is hard to explain within a simple physical picture or model, but surely it is closely related to the crystal surface because the most significant changes occur in the crystal surface. Two factors may play the major role: (1) The structure of the crystal surface. It promotes the ordered arrangement of ions, thus destroying the hydration structure and leads to the redistribution of water molecules. (2) The Coulomb field formed by anions and cations in the crystal makes the water molecules at the surface to rearrange with a specific orientation. These two factors, together with the hydration structure of ions, lead to the ordering distribution of water molecules in the FHS and slight disordering in dipole orientations.

As final remark, it should be point that, in the solution, the water molecules appear uniformly around the ions in the means of time average; however, the ordered arrangement of water molecules around an ion at the interface is more likely pinned. This is because the diffusion of ions and water molecules at the interface is much slower than that in the solution.

Summary

The solid–liquid interface between a NaCl crystal and its solution along the (100) direction of crystal has been studied systematically by using the molecular dynamics method. By calculating the various distributions of structural parameters, especially by distinguishing coordination environments of ions, we obtain the following main conclusions: (1) In the first layer nearest to the crystal surface, both cation and anion deviate from the epitaxial position of lattice. The magnitude of deviations strongly depends on the local coordination environments of ions. (2) The coordination between anions and cations has a little effect on the FHS of Na+ but has a significant impact on all aspects of the FHS of Cl–. (3) As approaching to the interface, the hydration structure of ions tends to be ordered; however, it is unconventional that the dipole orientation of water molecules in the FHS has the tendency of disordering. The current results may be helpful to further understand the thermodynamic and dynamic behavior of ions at interfaces.

Computational Details

In this work, the interfacial structure between NaCl crystal and its solution has been studied, in which the interface is parallel to the (001) surface of NaCl crystal. The SPC/E model[50] and the Joung–Cheatham (JC) model[51] were used for water molecules and ions, respectively. The JC model has been improved to reproduce the experimental hydration free energies of ions and the lattice constants.[51,52] The cross terms of the interaction potentials were calculated by the Lorentz–Berthelot combination rules. The detailed potential parameters are listed in Table .

Table 2

Parameters of the Joung–Cheatham Model and the SPC/E Model

	ε (kcal/mol)	σ (Å)	charge (e)
Na+	0.3526	2.166	+1
Cl–	0.0128	4.830	–1
O	0.1554	3.166	–0.8476
H	0	0	0.4238

All simulation in this work was performed by the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) software package,[53] which is a widely used molecular dynamics (MD) code. The MD time step was taken as 1 fs, and the temperature in all simulations was set as 298 K.

To set up the simulation system, the following steps were performed: (1) based on the saturated concentration (see below), a certain number of water molecule and the corresponding number of anions and cations were chosen. The ions were randomly distributed into water to form a solution. (2) At 298 K and 1 bar, the solution was equilibrated through long-time molecular dynamics simulation. (3) The solution and a NaCl crystal are combined to build the interface system. Then, an extensive simulation with the NPT ensemble was performed for 5 ns to equilibrate the system, in which the pressure was set as 1 bar with the damping coefficient being 1000 time steps. Finally, a further 25 ns simulation at 298 K with NVT ensemble was used to calculate the physical properties.

A snapshot of simulation system is shown in Figure , where three coordinate axes, , , and , are defined along the [100], [010], and [001] directions of NaCl crystals, respectively. The solid–liquid interface is parallel to the – plane, i.e., perpendicular to the direction. In the current work, the origin of axis is defined as the position of the outermost crystal layer. The system (57.8 Å × 57.8 Å × 151.2 Å) is composed of NaCl crystal (4000 Na+ and 4000 Cl–) and NaCl solution (1034 Na+, 1034 Cl–, and 9272 water molecules). The concentration of the solution is 6.2 mol/kg, which is the saturation concentration given by the JC model at room temperature (298 K)[17] and close to the experimental value (6.16 mol/kg at 298.15 K).[17] The saturation concentration has been thoroughly tested by others[54] and us. In fact, in our simulation of up to tens of nanoseconds, the system does equilibrate at this concentration of room temperature.

Figure 7

Snapshot of the interface system. The blue, magenta, green, and red spheres represent Na+, Cl–, H, and O, respectively. The , , and axes are defined along the [100], [010], and [001] directions of NaCl crystals, respectively. The solid–liquid interface is parallel to the plane, and perpendicular to the direction.

To calculate the number density distribution along the direction perpendicular to the interface ρ(z), the system was equally divided into small bins with the width of 0.2 Å along the direction. ρ(z) readswhere the angle brackets indicate the ensemble average, and N, z, and ΔV denote the number of atoms in the ith bin, the central position of this bin, and the volume of a bin, respectively. We have tested a few different widths of bin (from 0.2 to 0.6 Å) and have not observed any evident effect on final results.

Regarding the ionic hydration structure, the water in the first hydration shell (FHS) has the most important influence on the structure and properties of the system;[48,49] thus, our studies will only focus on FHS. The schematic illustration of FHS of an ion is shown in Figure . The cutoff radius (Rc) for FHS is defined as first minimum of radial distribution functions (RDFs) between ions and oxygens. We have found that Rc for Na+ and Cl– is 3.25 and 3.75 Å, respectively. To further characterize the local environment of an ion, two coordination numbers (nc and nw) are defined, where nw denotes the number of water molecules in FHS, the so-called hydration number, and nc is the number of the nearest counter ions. When calculating nc, the ions in the crystal are not taken into account.

As approaching to the interface, the spatial translational and rotational symmetries of ions in the solution may be broken partially. To characterize these effects, four angles (α, β, γ, and φ), which have been marked in Figure , are defined based on two unit vectors (1 and 2). 1 is the unit vector from the central ion to oxygen atoms, and 2 is the unit vector along the dipole of a water molecule. α, β, and γ were defined as the angle between 1 and three axes (, , and ), respectively. φ is the angle between 1 and 2. The distribution of α, β, γ, and φ was calculated as followswhere z and N are the positions of the ith bin and the number of ions in this bin, respectively. nθ(cos(θ)) denotes the number of water molecules, whose cos(α) (or cos β), cos(γ), and cos(φ) equal cos(θ). It needs to be pointed out that, due to the symmetry of the current system, the distribution of cos(β) is the same as that of cos(α).