Molecular dynamics study of salt–solution interface: Solubility and surface charge of salt in water Kazuya Kobayashi, Yunfeng Liang, Tetsuo Sakka, and Toshifumi Matsuoka Citation: The Journal of Chemical Physics 140, 144705 (2014); doi: 10.1063/1.4870417 View online: http://dx.doi.org/10.1063/1.4870417 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/14?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Solubility of NaCl in water by molecular simulation revisited J. Chem. Phys. 136, 244508 (2012); 10.1063/1.4728163 Molecular dynamics study of solubilization of immiscible solutes by a micelle: Free energy of transfer of alkanes from water to the micelle core by thermodynamic integration method J. Chem. Phys. 133, 074511 (2010); 10.1063/1.3469772 Molecular dynamics simulations of nonpolarizable inorganic salt solution interfaces: NaCl, NaBr, and NaI in transferable intermolecular potential 4-point with charge dependent polarizability (TIP4P-QDP) water J. Chem. Phys. 132, 024713 (2010); 10.1063/1.3269673 Solubility of KF and NaCl in water by molecular simulation J. Chem. Phys. 126, 014507 (2007); 10.1063/1.2397683 Solubility of KF in water by molecular dynamics using the Kirkwood integration method J. Chem. Phys. 117, 4947 (2002); 10.1063/1.1498820
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 80.167.90.33 On: Sun, 18 May 2014 11:52:52
THE JOURNAL OF CHEMICAL PHYSICS 140, 144705 (2014)
Molecular dynamics study of salt–solution interface: Solubility and surface charge of salt in water Kazuya Kobayashi,1 Yunfeng Liang,1,a) Tetsuo Sakka,2 and Toshifumi Matsuoka1,a) 1 2
Environment and Resource System Engineering, Kyoto University, Kyoto 615-8540, Japan Department of Energy and Hydrocarbon Chemistry, Kyoto University, Kyoto 615-8510, Japan
(Received 26 July 2013; accepted 24 March 2014; published online 14 April 2014) The NaCl salt–solution interface often serves as an example of an uncharged surface. However, recent laser-Doppler electrophoresis has shown some evidence that the NaCl crystal is positively charged in its saturated solution. Using molecular dynamics (MD) simulations, we have investigated the NaCl salt–solution interface system, and calculated the solubility of the salt using the direct method and free energy calculations, which are kinetic and thermodynamic approaches, respectively. The direct method calculation uses a salt–solution combined system. When the system is equilibrated, the concentration in the solution area is the solubility. In the free energy calculation, we separately calculate the chemical potential of NaCl in two systems, the solid and the solution, using thermodynamic integration with MD simulations. When the chemical potential of NaCl in the solution phase is equal to the chemical potential of the solid phase, the concentration of the solution system is the solubility. The advantage of using two different methods is that the computational methods can be mutually verified. We found that a relatively good estimate of the solubility of the system can be obtained through comparison of the two methods. Furthermore, we found using microsecond time-scale MD simulations that the positively charged NaCl surface was induced by a combination of a sodiumrich surface and the orientation of the interfacial water molecules. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4870417] I. INTRODUCTION
Aqueous solutions are fundamental in biology, chemistry, geology, and engineering.1–12 They have continuously attracted interest owing to the specific interactions between water and ions, such as ionic hydration,1 ion pairing,2 and hydrogen bond formation (or breaking).3 For sodium chloride, the surface (interface) in solution often serves as an example of an uncharged surface. However, laser-Doppler electrophoresis has shown some evidence that sodium chloride crystals have a positive charge in saturated solution.4 Experimental measurement of the NaCl–solution interfacial structure is still challenging because of the buried nature of the liquid–solid interface and the difficulty of controlling the surface morphology in solution. Experiments of the NaCl–solution interface, such as surface X-ray diffraction and atomic force microscope (AFM), are still limited to the NaCl surface with a few layers of water.13, 14 Molecular dynamics (MD) simulations have been performed to study the static and dynamic phenomena of the NaCl–solution interface. The mechanism of the surface charge of the salt is still under debate.15–26 On one hand, the surface charge is attributed to water having a specific orientation at the interface,15–19 which is based on the calculations with both classical15–17 and ab initio MD approaches.18, 19 On the other hand, it has been reported that sodium ions are more easily adsorbed on the sodium chloride crystal surface at the early stages of crystallization.20–23 a) Authors to whom correspondence should be addressed. Electronic ad-
dresses:
[email protected] and
[email protected]. kyoto-u.ac.jp
0021-9606/2014/140(14)/144705/8/$30.00
Furthermore, it has been reported that positively charged clusters have a longer lifetime than negatively charged or neutral clusters.24 In addition, research on dissolution has shown that chloride ions present a greater tendency to be hydrated than sodium ions.25, 26 Despite numerous simulation results, it is still not clear why the sodium chloride crystal is positively charged at the equilibrated stage, because all of the simulation results only reported the early stage of crystallization or dissolution, presumably owing to computational cost. Moreover, as mentioned above, a unified picture on the surface charges of the NaCl–solution interface has not been achieved. The solubility is defined as the point when the chemical potentials of a certain component are equal in solution and in the pure solid phase. Thus, we can calculate the solubility by calculating the chemical potential. Considerable effort has been made to calculate the chemical potential and determine the solubility using molecular simulations (both Monte Carlo and MD). Cui and Harris27 were the first to propose a method to calculate the solubility of sodium chloride with MD by estimating the chemical potentials of sodium chloride. Their basic idea was to separately calculate the chemical potentials in solution and in the solid phase and compare the values. Based on this idea, Ferrario et al. calculated the solubility of potassium fluoride using a special technique to treat the ionic charge.28 Later, solubility calculations from molecular simulations using thermodynamic approaches were improved by other researchers.29–31 Recently, Moucka et al.33, 34 have developed the osmotic ensemble Monte Carlo (OEMC) method to calculate solubility. Another way to calculate solubility is the direct method (kinetic approach).31, 35 Using the
140, 144705-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 80.167.90.33 On: Sun, 18 May 2014 11:52:52
144705-2
Kobayashi et al.
direct method, the interfacial system consisting of the pure crystal and its oversaturated solution is prepared. After equilibrating, the saturated solution is obtained, which gives the solubility value. This is the simplest and easiest way to determine the solubility. However, the reported solubility values of NaCl differ from 3.7 to 7.27 mol/kg even using same force field.31–35 It is not clear which approach gives the most reliable value of the solubility. Thus, validation of the different methods is required. In summary, molecular simulations (mainly MD) have been used in the past to investigate the interfacial structures of the NaCl–solution system. However, a fully equilibrated system has never been addressed in previous studies. In addition, a consensus about the value of solubility for a given force field has not yet been reached. In this study, for the solubility calculations, we performed two MD simulations: one is the interfacial system calculation, called the direct method calculation, and the other is the chemical potential calculation. The values of the solubility from the two approaches are compared. This will help us estimate the solubility of a given force field more confidently. For the fully equilibrated system, we performed microsecond time-scale molecular dynamics. We consider that the interfacial system is equilibrated if the concentration of the salt in solution is converged to a certain value within ∼1 μs. The interfacial structure of the equilibrated system and the surface charge of the salt surface in solution are then analyzed and discussed. II. COMPUTATIONAL METHODS
MD simulations were performed using GROMACS (version 4.5.4)36 for both the direct method and the free energy calculations. Using the same package for both calculations avoids any spurious effects induced by small differences in the computational protocols of different software. The SPC/E37, 38 model was used for water. Both the CHARMM27r force field39 and the Joung–Cheatham (JC) ionic model40 were used for the sodium and chloride ions. Our previous study shows that the calculated interfacial tension of SPC/E water and CHARMM27r hydrocarbons are in excellent agreement with the experimental data.41 It is therefore interesting to test whether the CHARMM27r force field is also accurate enough for the ions. The JC potential parameters have been optimized to reproduce the experimental hydration free energies of ions in water,42 ion–water energies, and the lattice constants and lattice energies of salt crystals.35, 40 It has been shown that this set of parameters reproduce well various additional experimental properties, such as activity coefficients, diffusion coefficients, residence time of atomic pairs, association constants, and solubilities.35 The potential parameters of the CHARMM27r force field and the JC model are summarized in Table I. It should be noted that the potential that we used is a two-body additive force field, for which the polarizable contribution is excluded. The particle mesh Ewald summation43 was used to treat the electrostatic interactions and the cutoff length for the van der Waals interaction was 10 Å. The standard long-range correction for dispersion was included. The Nose–Hoover thermostat44 and Parrinello–Rahman barostat45 were used for the NPT simu-
J. Chem. Phys. 140, 144705 (2014) TABLE I. Parameters of the Joung–Cheatham40 and CHARMM2739 force fields. Model
Ion
ε (kJ/mol)
σ (nm)
Charge (e− )
Joung and Cheatham
Na+ Cl− Na+ Cl−
1.476 0.0535 0.196 0.628
0.216 0.483 0.243 0.405
1.0 − 1.0 1.0 − 1.0
CHARMM27
lations. Anisotropic scaling without shear deformation was used for the direct method, that is, the box was always orthorhombic.
A. Direct method
The salt–solution interfacial system was used to calculate the solubility. First, we prepared the sodium chloride crystal and its 120% oversaturated solution from experimental results, which consisted of 500 sodium chloride pairs and 2000 water molecules with 320 sodium chloride pairs, respectively (Fig. 1). Aragones et al.31 performed MD simulations by placing pure NaCl in contact with pure water. However, this approach was not successful for solubility calculations as no single ion from the solid moved into the fluid phase even after 500 ns owing to the high free energy barrier for dissolution from a perfect crystal.26 The same phenomenon occurred for ab initio MD simulations for tens of picoseconds.18, 19 We suggest that a similar phenomenon might occur if the solution concentration is less than (even equal to) the solubility limit. Therefore, an oversaturated solution rather than a dilute solution system was used in the direct method calculation in this study. When the combined interface system was prepared, we started with a vacuum space of ∼1.0 nm between the crystal and the solution for the stability of the calculation. This vacuum space disappeared around 30 ps into the simulation. A microsecond-scale simulation was performed to equilibrate the systems. The system was judged to be equilibrated based
FIG. 1. Snapshots of the NaCl–water interface system at 298 K and 0.1 MPa taken at (a) the initial step and (b) after 4150.0 ns. The solution area is indicated by the transparent yellow region. The equilibrated box sizes for the x, y, and z directions were 2.89, 2.89, and 11.5 nm, respectively.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 80.167.90.33 On: Sun, 18 May 2014 11:52:52
144705-3
Kobayashi et al.
J. Chem. Phys. 140, 144705 (2014)
FIG. 2. (a) Density profile of NaCl–water interface system averaged for the final 1 ns of the simulation and (b) temporal evolution of concentration of ions in the solution area.
on the concentration of ions in the solution area and number densities from the interfaces (Fig. 2). After equilibrium, the solubility value was calculated from the number density of ions and water in the solution area, which was a 3 nm wide region in the z direction at the center of the solution, and the interfacial structure was analyzed from the final 1 ns trajectory, unless specified otherwise. The calculated solubility value (from the direct method) was then compared with the value calculated by the free energy method, which is described below. B. Chemical potential calculation
For the chemical potential calculation, we used the approach of Sanz and Vega30 and Aragones et al.31 for the solid and in-solution phases, respectively. The Einstein crystal46 and Lennard–Jones (L-J) fluid47 were used for the reference states in the solid phase and in-solution phase, accordingly. Note that, Sanz and Vega’s calculations were with the MC method and ours were with MD using the same code as the direct method (GROMACS 4.5.4).36 As mentioned in the Sec. I, we need to separately calculate the chemical potentials of the solid and in-solution phases. The free energy difference is calculated by thermodynamic integration (TI) using the coupling parameter λ, which controls the potentials of the two states. During TI, the potential energy of the system is described by following equation: Usystem = (1 − λ) UA + λUB .
(1)
The free energy difference between states A and B can be calculated by integrating the ensemble averages of the partial derivative of potential energy with respect to λ for a set of λ: 1 ∂Usystem dλ, (2) A = ∂λ 0 λ where the system behaves as state A when λ = 0. The system completely changes to state B when λ = 1. We equilibrated the system for 300 ps at a certain λ value and then performed a calculation for a further 300 ps at the same λ for analysis.
The final configuration of the previous calculation was used for the equilibration run at the next λ. The chemical potential of the solid phase was calculated based on the Gibbs free energy difference between the pure sodium chloride crystal and the Einstein crystal with TI. First, an ideal sodium chloride crystal with 500 NaCl pairs was prepared. An equilibration calculation with the NPT ensemble was performed for 0.5 ns to obtain the lattice parameters and average volume at ambient pressure and temperature. After equilibrating, the harmonic spring interaction was gradually switched on with the NVT ensemble using the average volume with the interatomic interaction. The interatomic interaction was then switched off with the NVT ensemble to obtain the pure Einstein crystal, for which the free energy value can be theoretically calculated. A total of 16 points (eight points for each step) were obtained. For each step, eight points were used for TI with the Gaussian quadrature.30 From this operation, the Helmholtz free energy difference between the pure ideal crystal and the pure Einstein crystal was calculated. The Gibbs free energy difference between the two crystals is described by Gsolid = Asolid + pVsolid , where Asolid is the free energy difference calculated by the previous operation and Vsolid is the average volume of the system from the NPT calculation. We thus obtained the Gibbs free energy difference between the Einstein crystal and the pure sodium chloride crystal. The free energy value of the Einstein crystal was calculated from the theoretical equation.30 From this, we calculated the free energy value of the pure sodium chloride crystal, Gsolid . The chemical potential of sodium chloride in the solid phase was then obtained from: μsolid = Gsolid /NNaCl , where NNaCl is the number of pairs of sodium and chloride. In this study, NNaCl = 500 was chosen for the chemical potential calculation of the solid phase. Sanz and Vega reported that there was no size or integration direction dependence in this system size regime.30 The chemical potential calculation of the in-solution phase was performed with TI between various concentrations of sodium chloride solution and L-J fluid, because the chemical potential of the in-solution phase is dependent on
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 80.167.90.33 On: Sun, 18 May 2014 11:52:52
144705-4
Kobayashi et al.
J. Chem. Phys. 140, 144705 (2014)
TABLE II. Solubility value of all simulations (experimental value: 6.16 mol/kg at 298.15 K49 ). SPC/E compatible ion System size This work
Large Small
Joung and Cheatham (2009)35 Aragones et al. (2012)31 Moucka et al. (2011)33 Moucka et al. (2013)34
Direct calc. (mol/kg)
Chemical potential calc. (mol/kg)
6.20 8.07 7.27 5.5 ... ...
2.79–6.52
4.8 4.8 3.7 CHARMM27
System size Small
This work
Direct calc. (mol/kg) 2.57
concentration. The system contained 500 water molecules and 5–70 pairs of sodium chloride, where the concentration range was 0.56–7.78 mol/kg. The operation was similar to the solid phase case: the system was first equilibrated with the NPT ensemble for 10 ns and then the TI was used with the NVT ensemble using the volume from the 2nd-order polynomial fitting of the average number density of sodium chloride of the previous equilibration run (0.56–7.78 mol/kg), except that 21 independent λ points were used instead of 16. We used the scheme proposed by Aragones et al. in 2012.31 That is, the free energy was divided into an ideal and a residual contribution. First, we obtained the derivative of the ideal part numerically using the number density of solutions at different concentrations. Second, we obtained the residual part by taking into consideration the free energy difference between the solution and the L-J fluid. The total chemical potential was the sum of the ideal and residual chemical potentials. This scheme is very useful to introduce curvature in the relationship between the chemical potential of NaCl and the concentration. To calculate the residual free energy, the free energy of the L-J fluid was chosen as the reference state, which can be calculated from the equation of state.47 The free energy difference was calculated by Gres = Asolution + pVsolution , where Asolution is obtained from the TI and Vsolution is the fitting result of the average volumes. The residual Gibbs free energy value of various concentrations can be obtained. The residual chemical potential of the in-solution phase is then described as μres =
∂Gres ∂NNaCl
.
(3)
Nwater ,p,T
Because we need to calculate the partial derivative of the Gibbs free energy with respect to the number of sodium chloride, we fitted the Gibbs free energy data by a polynomial function. Note that this method only works well when considering a narrow range of concentrations (see the supplementary material for further discussion).30, 48 Once the chemical potential of the solid phase and the chemical potential function of the in-solution phase were obtained, we determined the point when the solid and in-solution phases are equal, which is the solubility of sodium chloride.
Chemical potential calc. ...
III. RESULTS AND DISCUSSION A. Solubility
First, the solubility using the JC model calculated by the direct method is discussed. Figure 2 shows the number density profile of the salt–solution interface system and the temporal evolution of the concentration of ions in the solution area at 298 K and 1 bar. It took about 3.1 μs for the interface system to reach its equilibrium state (Fig. 2(b)). Simulation for a further microsecond was performed to confirm that the system was in equilibrium at least at the microsecond scale. The equilibrium concentrations of Na+ and Cl− in the system were almost identical, as expected. Thus, the average number density of each component in the solution area (over the last 100 ns) was used, giving a solubility of 6.20 ± 0.23 mol/kg. Table II shows that the calculated solubility value is located in the middle of the values reported by Aragones et al.31 (5.50 mol/kg) and Joung and Cheatham (7.27 mol/kg).35 Note that, all of the three calculations used the direct method and all the calculations were in the order of microseconds or half a microsecond. Because the simulation in this work is one of the longest calculations, we believe that 6.20 ± 0.23 mol/kg is close to the true solubility value of the JC ionic model in SPC/E water. It is remarkable that this value is in excellent agreement with experimental data (6.16 mol/kg at 298.15 K).49 Nevertheless, we caution that some long-time processes (longer than microseconds) might be missing in the direct method of MD simulations. Next, we calculated the chemical potentials of the solid and in-solution phases. In the solid phase case at 298 K, the harmonic springs were first introduced to the system of pure sodium chloride (Fig. 3(a)). The interatomic potentials (L-J parameters and atomic charge) were gradually switched off (Fig. 3(b)). From the two TIs, the Gibbs free energy difference between pure sodium chloride and the pure Einstein crystal was obtained. The calculated chemical potential for NaCl solid with the JC force field was −770.5 kJ/mol, which is in good agreement with the experimental value of around −770.8 kJ/mol50 and a previously calculated value with the MC method and the same force field (−770.9 kJ/mol).31 In the in-solution phase case, we used the L-J fluid as the reference state, which had L-J parameters of εref /kB = 78.2 K and σ ref = 3.14 Å. Simultaneously, the L-J parameters of all
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 80.167.90.33 On: Sun, 18 May 2014 11:52:52
144705-5
Kobayashi et al.
J. Chem. Phys. 140, 144705 (2014)
FIG. 4. Derivative of the Hamiltonian with respect to lambda (λ) for the transformation of the Joung–Cheatham NaCl solution into a Lennard–Jones fluid for solutions with various concentrations.
FIG. 3. (a) Derivative of the Hamiltonian as a function of lambda (λ) at 298 K calculated in the solid system: (a) the result of the first step (transition from real solid to Einstein-crystal with ionic interaction) and (b) the result of the second step (transition into pure Einstein-crystal). See text for definition of lambda (λ).
of the particles (except for hydrogen in water) were changed to those of the reference state (Fig. 4) and the charges of all of the particles were switched off during the TI. This calculation used a constant number of water molecules and various concentrations of sodium chloride. By implementing Simpson’s integration with Eq. (2), we obtained Asolution . Finally, we obtained the residual Gibbs free energy as function of the number of sodium chloride pairs in solution. This must be fitted to a polynomial function because the chemical potential of sodium chloride in solution is given by the partial derivative of the Gibbs free energy with respect to the number of sodium chloride pairs. After fitting a 2nd-order polynomial to the range of values from 5.77 to 7.78 mol/kg, the residual chemical potential as function of the number of sodium chloride (or concentration) was obtained as μres (NNaCl ) = 0.4976NNaCl − 767.19. The number density as function of the number of sodium chloride pairs as ρ(N/Å3 ) = 0.0334603 + 7.49909 × 10−5 × NNaCl − 4.2901 × 10−7 2 × NNaCl . The calculated solubility at 298 K was 6.42 mol/kg, which is very close to the experimental solubility (6.16 mol/ kg at 298.15 K)48 and our calculated value by the direct method (Table II). Figure 5 also shows the chemical poten-
tials of NaCl in solution after fitting a 2nd-order polynomial to the full range of data (i.e., 0.56–7.78 mol/kg). The calculated solubility is about 2.79 mol/kg, which is considerably less than the value calculated by the direct method. This indicates that the free energy calculations only work well when considering a narrow range of concentration. Furthermore, it seems that the dilute regime dominates the trend in the fitting process. As shown in Fig. S1 of the supplementary material,48 the curve is very similar if the fitting window consists of the dilute regime. Comparing with previous work (Fig. 5), the chemical potential of NaCl in solution is significantly different in dilute solution, but converges at high concentrations. In the OEMC method used by Lisal et al.,32 the chemical potential is an input parameter and the concentration of the ions is an output of the simulation. OEMC is an open ensemble simulation and has advantages in the calculation of the chemical potential in the low concentration
FIG. 5. Chemical potentials of the in-solution phase at 298 K as a function of sodium chloride concentration. For comparison, the chemical potential calculated by Aragones et al.31 and Moucka et al.33 are also shown. Note: −386.8 kJ/mol is added to the calculated chemical potentials. The justification for the shift is discussed in Ref. 31. The residual chemical potential is assumed to be a linear function of the salt concentration in the solution. This assumption only holds for a narrow range of concentrations.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 80.167.90.33 On: Sun, 18 May 2014 11:52:52
144705-6
Kobayashi et al.
range.32 Similar to Aragones et al.’s MC method,31 our MD simulation uses a closed ensemble and has disadvantages in the low concentration range. Nevertheless, we suggest a solubility value of ∼6.42 mol/kg by comparison with the direct method calculations. If instead we rely on the chemical potential calculations only (not biased by the direct method calculation), we suggest a solubility value from 2.79 mol/kg to 6.52 mol/kg. It is shown in the supplementary material48 that the solubility variation is partially due to the fact that the plateau of the chemical potential in the in-solution phase is very close to that in solid phase, and partially due to the statistical uncertainty of the calculated data (i.e., the calculated residual Gibbs free energy). For the direct method, various additional systems were used to investigate the sensitivity of the solubility. In the first case, the system had 256 sodium chloride pairs for the crystal and its 120% oversaturated solution had 1208 water molecules and 190 sodium chloride pairs. The calculated concentration value was ∼8.07 mol/kg, which is much larger than the values calculated by the direct method and the free energy calculations. Because of the large difference, we suggest that this system (with 256 sodium chloride pairs for the crystal) is too small for the solubility calculation. An interface system with 500 sodium chloride pairs for the crystal but with surface defects (voids) gave a calculated solubility value of 7.08 mol/kg, which is slightly larger than the values mentioned above but close to the previously calculated value of 7.27 mol/kg by Joung and Cheatham.35 Therefore, these small differences likely result from the difference of the initial configurations. In fact, the relaxation time of the latter system (with some voids on the salt surface) was about 500 ns, which is very close to relaxation times in previous calculations.35 We also used the CHARMM27r force filed for the ions. The size of the system was similar to the small system mentioned above (i.e., with 256 sodium chloride pairs for the crystal). The calculated results (∼2.57 and 2.97 mol/kg from the free energy and direct methods, respectively) are very different from the experimental value. Given that the small system tends to overestimate the solubility limit in the direct method, we suggest that the CHARMM27r ion models should not be used in calculations of interfacial systems. In fact, a preliminary result with the CO2 /brine/silica system has shown that salt crystal forms on the silica surface far below the experimental solubility limit.
B. Interfacial structure
The solubility with the JC model showed better agreement with the experimental data. Because the system was assumed to be equilibrated, we performed further analysis of the interfacial structure with the JC model for the large system (the system without voids on the salt surface). The crystal surface step structure at the end of simulation (for the final 1 ns) is shown in Fig. 6. The crystal surface was defined based on the distance between the initial crystal surface and ions in solution adjacent to the surface. If the average distance during the final 1 ns of the simulation was less than the first minimum (0.36 nm) of the Na–Cl radial distribution function (RDF) in
J. Chem. Phys. 140, 144705 (2014)
FIG. 6. Surface step structure for the final 1 ns of the simulation. Note: NLayer = 0 means no crystal layer formed from the initial crystal.
the solution system, the ion was defined as a new crystal surface and the crystal surface was renewed. The crystal surface determination was performed for all ions in proximity to the salt surface. Figures 7(a) and 7(b) show the temporal evolution of Na+ and Cl− ion number densities for 0.25 nm blocks at different distances from the left and right surfaces, respectively. It is important to confirm that the direct method was indeed equilibrated. As shown in Fig. 7(a), the first newlyformed crystalline block (0.25–0.5 nm) of the left surface was equilibrated at around 1.2 μs, while the first newly-formed crystalline block (0.25–0.5 nm) of the right surface (Fig. 7(b)) was equilibrated at around 3.2 μs. Furthermore, all of the blocks (0.25–1.0 nm) were equilibrated at around 3.4 μs. From Figures 7(a) and 7(b), the number density of sodium ions at both surfaces (i.e., the edge of the solid phase) is always higher than that of the chloride ions. Initially, the number density of sodium ions at a distance of 0.5 nm is higher than that of the chloride ions. After the layer forms at a distance of 0.5 nm, the number density of sodium ions at a distance of 0.75 nm becomes higher than that of chloride ions. Therefore, it is suggested that the sodium chloride crystal surface is positively charged, which agrees with the experimental findings.4 In contrast, the number density of chloride ions at a distance of 1.0 nm is higher than that of sodium ions. This can be regarded as counter-ion behavior in proximity to a charged surface. This hypothesis can be clearly verified by Fig. 8 (averaged for the last 1 ns), in which the first peak intensity (from left) of the sodium ions at the interface is much higher than that of the chloride ions. In addition, we also investigated the solvent effect at the crystal surface by calculating the water orientation at the surface. The water orientation at the surface can be defined by two angles: the angle between the surface normal direction (z-vector) and bisector vector in water molecule, and the angle between the z-vector and oxygen–hydrogen vector. Figure 8 shows the results of the number density of the water molecule having a certain angle with the z-vector. In total, three water layers were found at the surface, from which we can construct the orientation of water molecules. The water molecules have one hydrogen atom pointing away from the crystal surface, as shown schematically in Fig. 8(c). Thus, we can assume that the water layer also makes the surface positively charged. The fact that chloride ion peak position (∼1.9 Å away from the first peak) is different from the crystal periodicity (∼2.9 Å) verifies the
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 80.167.90.33 On: Sun, 18 May 2014 11:52:52
144705-7
Kobayashi et al.
J. Chem. Phys. 140, 144705 (2014)
FIG. 7. Temporal evolution of number densities of the ions, where +X nm means the area between the distances X−0.25 and X nm from the initial crystal layer: (a) calculated from the front surface and (b) calculated from the other surface owing to the periodic boundary condition.
FIG. 8. Two-dimensional orientation profiles of water with respect to the NaCl surface normal and z-coordinate. The dipole moment of water (a) and OH bond of water (b) are selected as the molecule axis. The Na+ and Cl− number density profile (from Fig. 2(a)) is overlaid with the orientation profile for comparison. A schematic diagram of the interface structure (c) is drawn based on Fig. 6 and the two panels (a) and (b).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 80.167.90.33 On: Sun, 18 May 2014 11:52:52
144705-8
Kobayashi et al.
assumption (Fig. 8). That is, the solvent has a significant effect on the ionic structures at the interface. In summary, the positively charged sodium chloride surface consists of both sodium-rich steps at the surface and a positively charged water layer. This finding also holds for the surface starting with some voids, as mentioned in Sec. III. The only difference is that the time scale required to reach the equilibrium is shorter. IV. CONCLUSIONS
Solubility calculations with kinetic and thermodynamic approaches were in good agreement if the system size of the salt–solution interface system was sufficiently large and the simulation time was sufficiently long. The equilibrium time for the interface system takes from 500 ns to 3.4 μs depending on the initial configuration. Based on our simulation results, it can be assumed that both analyses can correctly predict the solubility with the ion and water models. However, care must be taken for the free energy calculations. In particular, choosing the fitting window of the polynomial fit for the in-solution phase. This can be made easier by comparing the two different approaches. Furthermore, the solubility can be incorrect with certain ion and water models, so the models need to be carefully chosen. Since the JC and SPC/E water models have a good solubility value and crystallization can be observed, the interfacial structure between the newly formed crystal and its saturated solution was analyzed. The concentration of sodium ions at the interface was high for almost all of the period of crystallization. We concluded that the positively charged surface of the sodium chloride crystal (in aqueous solution) is induced by a sodium-rich step structure and the water layer in proximity to the interface. ACKNOWLEDGMENTS
The authors acknowledge the financial support of the Japanese Society for the Promotion of Science (JSPS) through a Grant-in-Aid for Scientific Research A (No. 24246148), JOGMEC, JST/JICA-SATREPS, and JAPEX. We also wish to thank William R. Smith, Caetano R. Miranda, and Yasuhiro Fukunaka for valuable discussions. 1 H.
Ohtaki and T. Radnai, Chem. Rev. 93, 1157 (1993). Marcus and G. Hefter, Chem. Rev. 106, 4585 (2006). 3 Y. Marcus, Chem. Rev. 109, 1346 (2009). 4 J. D. Miller, M. R. Yalamanchili, and J. J. Kellar, Langmuir 8, 1464 (1992). 5 H. O. Yildiz and N. R. Morrow, J. Pet. Sci. Eng. 14, 159 (1996). 6 P. Zhang, M. T. Tweheyo, and T. Austad, Colloids Surf., A 301, 199 (2007). 7 S. Fathi, T. Austad, and S. Strand, Energy Fuels 24, 2514 (2010). 8 V. A. Tabrizy, A. A. Hamouda, and R. Denoyel, Energy Fuels 25, 1667 (2011). 9 V. A. Tabrizy, R. Denoyel, and A. A. Hamouda, Colloids Surf., A 384, 98 (2011). 10 T. Hassenkam, C. S. Pedersen, K. Dalby, T. Austad, and S. L. S. Stipp, Colloids Surf., A 390, 179 (2011). 2 Y.
J. Chem. Phys. 140, 144705 (2014) 11 T.
Hassenkam, A. C. Mitchell, C. S. Pedersen, L. L. Skovbjerg, N. Bovet, and S. L. S. Stipp, Colloids Surf., A 403, 79 (2012). 12 M. Hovland, H. Rueslatten, H. Johnsen, B. Kvamme, and T. Kuznetsova, Mar. Pet. Geol. 23, 855 (2006). 13 J. Arsic, D. M. Kaminski, N. Radenovic, P. Poodt, W. S. Graswinckel, H. M. Cuppen, and E. Vlieg, J. Chem. Phys. 120, 9720 (2004). 14 A. Verdaguer, G. M. Sacha, M. Luna, D. F. Ogletree, and M. Salmeron, J. Chem. Phys. 123, 124703 (2005). 15 H. Shinto, T. Sakakibara, and K. Higashitani, J. Phys. Chem. B 102, 1974 (1998). 16 H. Du and J. D. Miller, J. Phys. Chem. C 111, 10013 (2007). 17 S. Yamanaka, A. Shimosaka, Y. Shirakawa, and J. Hidaka, J. Nanopart. Res. 12, 831 (2010). 18 L. Liu, M. Krack, and A. Michaelides, J. Am. Chem. Soc. 130, 8572 (2008). 19 L. Liu, M. Krack, and A. Michaelides, J. Chem. Phys. 130, 234702 (2009). 20 H. Shinto, T. Sakakibara, and K. Higashitani, J. Chem. Eng. Jpn. 31, 771 (1998). 21 E. Oyen and R. Hentschke, Langmuir 18, 547 (2002). 22 Y. Yang and S. Meng, J. Chem. Phys. 126, 044708 (2007). 23 K. Kadota, Y. Shirakawa, M. Wada, A. Shimosaka, and J. Hidaka, J. Mol. Liq. 166, 31 (2012). 24 L. Degrève and F. da Silva, J. Chem. Phys. 111, 5150 (1999). 25 H. Ohtaki, N. Fukushima, E. Hayakawa, and I. Okada, Pure Appl. Chem. 60, 1321 (1988). 26 L. Liu, A. Laio, and A. Michaelides, Phys. Chem. Chem. Phys. 13, 13162 (2011). 27 S. Cui and J. Harris, J. Phys. Chem. 99, 2900 (1995). 28 M. Ferrario, G. Ciccotti, E. Spohr, T. Cartailler, and P. Turq, J. Chem. Phys. 117, 4947 (2002). 29 A. S. Paluch, S. Jayaraman, J. K. Shah, and E. J. Maginn, J. Chem. Phys. 133, 124504 (2010). 30 E. Sanz and C. Vega, J. Chem. Phys. 126, 014507 (2007). 31 J. L. Aragones, E. Sanz, and C. Vega, J. Chem. Phys. 136, 244508 (2012). 32 M. Lisal, W. R. Smith, and J. Kolafa, J. Phys. Chem. B 109, 12956 (2005). 33 F. Moucka, M. Lisal, J. Skvor, J. Jirsak, I. Nezbeda, and W. R. Smith, J. Phys. Chem. B 115, 7849 (2011). 34 F. Moucka, I. Nezbeda, and W. R. Smith, J. Chem. Phys. 138, 154102 (2013). 35 I. S. Joung and T. E. Cheatham III, J. Phys. Chem. B 113, 13279 (2009). 36 B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008). 37 H. J. C. Berendsen, J. Grigera, and T. Straatsma, J. Phys. Chem. 91, 6269 (1987). 38 J. Alejandre, D. Tildesley, and G. Chapela, J. Chem. Phys. 102, 4574 (1995). 39 J. Klauda, B. Brooks, A. MacKerell, Jr., R. Venable, and R. Pastor, J. Phys. Chem. B 109, 5300 (2005). 40 I. S. Joung and T. E. Cheatham III, J. Phys. Chem. B 112, 9020 (2008). 41 M. Kunieda, K. Nakaoka, Y. Liang, C. R. Miranda, A. Ueda, S. Takahashi, H. Okabe, and T. Matsuoka, J. Am. Chem. Soc. 132, 18281 (2010). 42 R. Schmid, A. M. Miah, and V. N. Sapunov, Phys. Chem. Chem. Phys. 2, 97 (2000). 43 T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). 44 S. Nose, Mol. Phys. 52, 255 (1984). 45 M. Parrinello and A. Rahman, J. Chem. Phys. 76, 2662 (1982). 46 D. Frenkel and A. J. C. Ladd, J. Chem. Phys. 81, 3188 (1984). 47 J. Kolafa and I. Nezbeda, Fluid Phase Equilib. 100, 1 (1994). 48 See supplementary material at http://dx.doi.org/10.1063/1.4870417 for details. 49 D. R. Lide, Handbook of Chemistry and Physics, Internet Version, 87th ed. (Taylor and Francis, Boca Raton, FL, 2007). 50 M. W. Chase, Jr., “NIST-JANAF thermochemical tables,” J. Phys. Chem. Ref. Data Monogr. No. 9 (American Institute Physics, Woodbury, NY, 1998).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 80.167.90.33 On: Sun, 18 May 2014 11:52:52