Dynamic and Static Properties of Aqueous NaCl Solutions at 25°C as a Function of NaCl Concentration: A Molecular Dynamics Simulation Study

Abstract

We present the result of molecular dynamics (MD) simulations to calculate the molar conductivity of NaCl in SPC/E water at 25^°C as a function of NaCl concentration (c) using Ewald sums employing a velocity Verlet algorithm. It is found that the MD result for Λ_m with Ewald sum parameter κ = 0.10 Å⁻¹ gives the closest one to the experimental data and that the obtained radial distribution functions (r) with κ = 0.10 Å⁻¹ show a dramatic change with a very deep minimum of (r) and, as a result, sharp maxima of (r) and (r) at the distance 9.95 Å, which indicates a characteristic of ionic atmosphere, the basis of the Debye–Hückel theory of ionic solutions. The static and dynamic properties of NaCl (aq) solutions are analyzed in terms of radial distribution functions, hydration numbers, coordination numbers around Na⁺ and Cl⁻, residence times of water around Na⁺ and Cl⁻, water diffusion, and ion-ion electrostatic energies to explain the behavior of the molar conductivity Λ_m of NaCl obtained from our MD simulations.

1. Introduction

The theory that an electrolyte is dissolved in water and dissociated into ions was proposed by Arrhenius (Arrhenius’s theory of electrolytic dissociation), who won the Nobel Prize in Chemistry in 1903, in his doctoral degree in 1884. Before this theory was raised, it was well known that a solution dissolved in an electrolyte passes electricity, and Faraday explained that the electrolyte dissolved in water is decomposed into charged ions by the action of this electricity. By measuring the molar electrical conductivity according to the concentration of various salt solutions, Kohlrausch obtained the relation (equation (1)) between the molar conductivity and the molar concentration (c) of a strong electrolyte (Kohlrausch’s law), and he found that each ion moved independently in an electric field. In the case of an electrolyte dissolved in a solvent and dissociated into cations and anions, the limiting molar conductivity in an infinitely dilute solution is related to the limiting molar conductivity ( and ) of each ion, as shown in equation (2) (Kohlrausch’s law of the independent migration of ions).where K is a constant.where the limiting molar conductivities of Na⁺ and Cl⁻ at 25^oC are measured as 5.01 and 7.63 mS m² mol⁻¹ experimentally [1].

In 1925, Onsager extended the Debye–Huckel theory of electrolyte solutions, published two years ago, to come up with the Debye–Huckel–Onsager theory, which became an integral part of his Nobel Prize winning achievements. Around some ions, ions with opposite charges are surrounded by electrostatic forces. In the absence of an external electric field, the distribution of counter ions around the central ion is symmetric. However, when the electric field acts, the ions at the center and the ions around them receive forces in opposite directions to create an asymmetric ionic atmosphere, and thereby, an electrostatic force acts in the direction opposite to the direction in which the ions move. The higher the concentration of the electrolyte, the greater this effect becomes, and thus, the molar conductivity becomes less than in a dilute solution. This reduction of the ions’ molar conductivity is called the relaxation effect.

Another effect of the ionic atmosphere is the electrophoretic effect. Ions are hydrated with solvent molecules, and these solvent molecules move with the ions. Since the central ions and the surrounding ions move in opposite directions, the central ions move against the flow of the solvent moving together with the surrounding ions. Onsager took these two points into consideration and treated the ion motion as a zigzag Brownian motion rather than a straight line in the electric field and derived the following molar conductivity equation. The Debye–Huckel–Onsager theory is an attempt to obtain quantitative expressions which leads to a Kohrrausch-like expression:where A and B are coefficients determined by temperature and the dielectric constant and viscosity of the solvent. The numerical values of A and B in water are 6.02 mS m² mol⁻¹/(mol·L⁻¹)^1/2 and 0.229 (mol L⁻¹)^1/2 [1]

A systematic conductivity study of sodium chloride-water + 1, 4-dioxane in dilute solutions (c < 0.01 mol-dm⁻³) [2] has been conducted covering a broad solvent composition range at temperatures from 5 to 35^°C. Accurate viscosity and permittivity data were determined for the organic solvent system. Evaluation of the limiting molar conductivity , ionic conductivities and, and the association constant is based on the chemical model [3] of electrolyte solutions, including short-range forces. In the next paper [4], they reported on an extension of this study to concentrated solutions, covering the temperature range from 5 to 35^°C at a solvent mole fraction of dioxane from xD = 0 to 0.65 and electrolyte concentration up to the limit of salt solubility in the respective mixtures and up to 5 mol-dm⁻³ in water. Data analysis is based on the mean spherical approximation (MSA) [5]. Comparison is made with the data representation by the empirical Casteel–Amis equation [6]. The association constants of the MSA are compared with those from chemical model calculations at low concentrations (lcCM). Later on, an investigation of ion-pairing of alkali metal halides in aqueous solutions using the electrical conductivity and the Monte Carlo computer simulation methods was reported. [7].

In our previous paper [8], we reported molecular dynamics (MD) simulations for the ionic mobilities μ of Na⁺ and Cl⁻ ions, in the system of Na⁺ (or Cl⁻ ) ion + 1023 water molecules, at 25°C using reaction field correction (RFC), simple truncation (ST), and Ewald sums employing Gear’s fifth order predictor-corrector and velocity Verlet algorithms. We found that MD simulations using Ewald sum with κ = 0.20–0.30 Å⁻¹ employing a velocity Verlet algorithm gives the best result for and at 25 C, even though the calculated overestimates the experimental data and the underestimates them. This behavior of is opposite to that reported in our other MD simulations [9].

In this study, we have carried out MD simulations to calculate the molar conductivity Λ_m of NaCl in SPC/E water at 25°C as a function of NaCl concentration (c) using Ewald sums employing a velocity Verlet algorithm. The primary goal of this study is to analyze radial distribution functions, hydration numbers, coordination numbers, residence times, and ion-ion electrostatic energies water diffusion to understand the behavior of the ionic atmosphere around Na⁺ and Cl⁻ ions through our MD simulations. This paper is organized as follows: Section 2 contains a brief description of molecular models and MD simulation methods followed by Section 3, which presents the results of our simulations. Our conclusion is summarized in Section 4.

2. Molecular Models and Molecular Dynamics Simulation Details

In order to establish systems of electrolyte solution of NaCl for MD simulation, we consider a system of c = 1 molar concentration (mol/L) of NaCl solution with 36 of Na⁺ and Cl⁻ ions at 25^°C. Then, the volume of this solution is simply equal to V = V_soln = 36000 cm³/N_A with N_A, Avogadro’s number, which is the fixed volume for all the current MD simulations of NaCl in electrolyte systems, which gives the length of simulation box L = 39.1009 Å. It is very hard to find the relation between the molar concentration (c) and the density of the solution (ρ_t, g/cm³) of NaCl solution in the literature. Instead, we use the experimental data of the weight percent (%W_NaCl) and the density (ρ_t, g/cm³) of NaCl solutions (Table 1) to find the numbers of ions (N_NaCl) and water (N_w). First, the relationship between ρ_t and %W_NaCl is obtained by using the least squares method from Table 1.

Table 1 

Experimental data between %W_NaCl and ρ_t at 25^oC. The relationship (ρ_t = 0.00758333853 %W_NaCl + 0.9942949) is obtained by using the least squares method.

For c = 1 mol/L of N_NaCl = 36 with V = 36000 cm³/N_A,

Also,

Substituting equations (5) and (6) into equation (4), we obtain N_w = 1956. Numbers of water molecules for given c are obtained using the same process and listed in Table 2.

Table 2 

The MD simulated systems in this study.

The SPC/E model [10] was adopted for water-water and ion-water. The pair potential between water and ion has the TIPS form [11]:where and are Lennard–Jones (LJ) parameters between oxygen on a water molecule and an ion i, q_j is the charge at site j in water, and q_i is the charge on ion i. Moreover, and are the distances between ion i and an oxygen site of a water molecule and between ion i and a charge site j in water, respectively. The LJ parameters σ_io and ε_io between oxygen on an SPC/E water molecule and the Na⁺ or Cl⁻ ion i were fitted to the binding energies of small clusters of ions by Dang et al. (σ_OO = 2.876 Å, σ_NaO = 2.876 Å, σ_ClO = 3.785 Å, ε_OO = 0.6502 kJ/mol, and ε_NaO = ε_ClO = 0.5216 kJ/mol) [12, 13]. The electrostatic charges of SPC/E water and ions are , , , and . The pair potential energy of ion-ion also has the TIPS form [11]:where and are Lennard–Jones (LJ) parameters between ion i and ion j, q_i is the charge on ion i, and is the distance between ion i and ion j. The LJ parameters σ_ij and ε_ij between ion i and ion j are given by σ_NaNa = 2.583 Å, σ_ClCl = 4.401 Å, σ_NaCl = 3.492 Å, and ε_NaNa = ε_ClCl = ε_NaCl = 0.4184 kJ/mol [12, 13]. Recently, even though many new models including Lennard–Jones parameters and the electrostatic charges of ions for aqueous electrolyte solutions have been developed [14–16], we have used the SPC/E model for water-water and ion-water, equation (7), since this model was used in our previous studies with and without Ewald sums [8, 12, 17].

Each MD simulation was carried out in the canonical ensemble (NVT fixed), and the density of NaCl solution is given in Table 1, which corresponds to a cubic box length of L = 39.1009 Å for all the MD simulations. The usual periodic boundary condition in the x-, y-, and z-direction and the minimum image convention for pair potential were applied. Gaussian kinetics [18, 19] was used to control the temperature, and a quaternion formulation [20, 21] was employed to solve the equations of rotational motion about the center of mass of rigid SPC/E water molecules. We have employed a velocity Verlet algorithm [22] for the time-integration algorithms: with a time step of 1 femtosecond.

For the Ewald sums, the Ewald sum parameters κ are varied from 0.05 to 0.25 Å⁻¹ in 0.05 Å⁻¹ increments, and the cutoff sphere of radius R_c in the real space and k_max in the reciprocal space are chosen as L/2 and 7, but |k|² ≤ 27. MD runs of 500,000 time steps each were needed for the ion-water system to reach equilibrium. The equilibrium properties were then averaged over 10 blocks of 100,000 time steps (0.1 ns) for a total of 1,500,000 (1.5 ns) for the systems of given 5 values of κ and given 9 molar concentrations (mol/L). The configurations of water molecules and Na⁺ (Cl⁻) ion(s) were stored every 5 time steps for further analysis.

There are two routes to calculate self-diffusion coefficients of water from MD simulations: the Einstein relation from the mean square displacement (MSD),and the Green–Kubo (GK) relation from the velocity autocorrelation (VAC) function,

3. Results and Discussion

3.1. Diffusion Coefficient

and at 25°C in SPC/E water were obtained from mean square displacements (MSDs), equation (9), of Na⁺ and Cl⁻ from our MD simulations using Ewald sum with κ = 0.05–0.25 Å⁻¹ employing a velocity Verlet algorithm. We have listed sum of Na⁺ and Cl⁻ diffusion coefficients in Table 3 and compared D as a function of molar concentration (c) of NaCl with the experimental data [23] in Figure 1. Figure 1 shows that the MD results D for all the Ewald sum parameter κ overestimate the experimental data (○) except those (◇) with κ = 0.10 Å⁻¹ at 3.0 and 4.0 mol/L and those (□) with κ = 0.05 Å⁻¹ at 2.0, 3.0, and 4.0 mol/L.

Table 3 

Diffusion coefficients (D, 10⁻⁹ m²·sec⁻¹) of NaCl in SPC/E water at 25°C employing a velocity Verlet algorithm, where n is the number of (Na⁺ or Cl⁻) ion(s).

Figure 1 

Comparison of the diffusion coefficients (D =  + ) of NaCl at 25°C as a function of the molar concentration (c) of NaCl for given Ewald sum parameters (κ) employing a velocity Verlet algorithm obtained from MSDs of our MD simulations with the experimental data [23].

When compared with the experimental data, the best MD results for D are obtained at 0.0278–0.50 mol/L with κ = 0.25 Å⁻¹ (▽), at 1.0 M with κ = 0.05 Å⁻¹ (□), and at 2.0–4.0 mol/L with κ = 0.10 Å⁻¹ (◇). This best result set of D is also listed in Table 3. In our previous MD study [8], the Ewald sum parameter κ = 0.20–0.30 Å⁻¹ gives the best result for the ionic mobility µ_Na+ and µ_Cl- at 25^°C in single-ion systems in 1023 SPC/E water molecules. However, in the current work for multi-ion systems, the numbers of SPC/E water molecules are in the range of N_w = 1800–2000. Note that the best results of D for N_NaCl = 1, 2, 4, 8, and 18 (the A1 ∼ A5 MD systems) are obtained with κ = 0.25 Å⁻¹ which also gives the best result for µ_Na+ and µ_Cl- in single-ion systems in 1023 SPC/E water molecules. [2].

In Figure 2, we have compared the best obtained D (□) and D with κ = 0.10 Å⁻¹ (◇) as a function of c with the experimental data (○). However, the obtained structural and dynamic properties of the MD systems are totally different for the different values of the Ewald sum parameter κ, which is discussed in the following sections. Instead of choosing the MD systems for the best obtained D to the experimental data with different κ, we have chosen those of the same κ = 0.10 Å⁻¹ which give D (◇) closest to the experimental data among those with the same κ, but less closer to the best obtained D with different κ. Figure 2 also shows and with κ = 0.10 Å⁻¹ as a function of c. D_Cl- at all the c (mol/L) is greater than D_Na+, and these two values become equal as c increases.

Figure 2 

Comparison of the best diffusion coefficients (D), those with κ = 0.10 Å-1, and and at 25°C as a function of c with the experimental data [23]. The inset shows comparison of molar conductivity (Λm) of NaCl as a function of c1/2 obtained from D with the experimental data.

At very low concentrations, the D values of NaCl are much higher than the experimental data, the values at 0.222 and 0.5 mol/L are slightly higher than the experimental data, and the values at 1.0 and 2.0 mol/L are similar to the experimental data. At 3.0 and 4.0 mol/L, it is slightly lower than the experimental data. This result can be highly dependent on the choice of the Ewald sum parameter, κ, as discussed above.

3.2. Molar Conductivity

We plotted molecular conductivity of NaCl as a function of c^1/2 obtained from D in our MD simulations in the inset of Figure 2, where , where is the charge on the ion in units of the electronic charge e, F is the Faraday constant, R is the gas constant, and T is the absolute temperature.

Both the experimental data and the MD results for of NaCl show two linear dependences (Ex1 and Ex2) as a function of c^1/2. The value of the first slope for the first 4 experimental points (Exp. Ex1) and that of the second slope for the rest 6 points (Exp. Ex2) are −8.83 and −2.94 mS m² mol⁻¹/(mol·L⁻¹)^1/2, respectively. Those for the first 3 MD results (MD Ex1) and the rest 6 results (MD Ex2) are −34.73 and −4.47 mS m² mol⁻¹/(mol·L⁻¹)^1/2, respectively. According to the Debye–Huckel–Onsager (DHO) theory, equation (3), the slope is equal to -8.91 mS m² mol⁻¹/(mol·L⁻¹)^1/2 with the numerical values of A and B in water, 6.02 mS m² mol⁻¹/(mol·L⁻¹)^1/2 and 0.229 (mol L⁻¹)^1/2, and of NaCl, 12.64 mS m² mol⁻¹ [1]. In conclusion, the MD results for of NaCl show two linear dependences like the experimental result, but the value of the first slope is too large for the DHO theory and the value of the second slope is similar to the experimental result.

3.3. Radial Distribution Functions and Running Hydration Numbers n_iw (r)

The radial distribution functions, (r) and g_ih (r), and the running hydration numbers, n_io(r), for the ions (Na⁺ or Cl⁻) and the O or H atoms of SPC/E water molecules are shown in Figure 3 for the A6 MD system (c = 1.0 mol/L, each Na⁺(Cl⁻) 36 ions with 1956 water molecules) with κ = 0.10 Å⁻¹ as a function of distance. The running hydration number is defined as

Figure 3 

Radial distribution functions _io(r) and _ih(r) and the running hydration numbers n_io (r) and n_ih (r) of SPC/E water molecules for the A6 MD system at 25°C as a function of the distance r (Å) between the ion (i) and the oxygen (O) and hydrogen atom (H) of a water molecule using Ewald sum with κ = 0.10 Å − 1.

Figure 3 shows the typical (r) functions with positions of _ClO (r) at the first maximum 3.25 Å and the first minimum 3.95 Å, of _Cl-H (r) at the first maximum 2.30 Å and the first minimum 3.10 Å, _NaO (r) at the first maximum 2.35 Å and the first minimum 3.15 Å and _NaH (r) at the first maximum 3.00 Å and the first minimum 3.75 Å. The positions of the first maxima and minima of _io (r) and g_ih (r) are almost the same for all the values of c and κ considered in this study and the magnitudes of those maxima and minima are very slightly different for all c and κ. These radial distribution functions _iw (r) of multi-ion systems are almost the same to those of single-ion systems [2].

The hydration number n_io (R₁) in the first shell is defined through equation.(11) from the ion-oxygen distribution functions _io (r) where the upper limit of integration R₁ is the radius of the first hydration sphere, which corresponds to the first minimum in _io (r). That is, R₁ = 3.95 Å for _ClO (r) and 3.15 Å for _NaO (r). We have listed the average hydration numbers n_NaO (R₁) and n_ClO (R₁) for the values of κ = 0.05, 0.10 and 0.25 Å⁻¹ in the first hydration shells of Na⁺ and Cl⁻ in SPC/E water at 25°C in Table 4. The values of n_NaO (R₁) and n_ClO (R₁) are almost the same for κ = 0.05 Å⁻¹ (5.64–5.85 and 7.55–7.62), κ = 0.10 Å⁻¹ (5.67–5.90 and 7.65–7.82), and κ = 0.25 Å⁻¹ (5.47–5.69 and 7.51–7.68).

Table 4 

The average hydration numbers (n_io (R₁)) and residence times (τ_io, ps) of water molecules in the first hydration shells of Na⁺ and Cl⁻ in SPC/E water at 25°C employing a velocity Verlet algorithm where R₁ = 3.15 Å for Na⁺ and R₁ = 3.95 Å for Cl⁻.

3.4. Radial Distribution Functions _ii (r) and Running Coordination Numbers n_ii (r)

The MD systems of small numbers of ion(s) – the A1 [each Na⁺(Cl⁻) 1 ion] and the A2 [each Na⁺(Cl⁻) 2 ions] show unstructured _ii(r) functions due to the very small numbers of ion(s). The radial distribution functions, _NaNa (r), _NaCl (r), and _ClCl (r), and the running coordination numbers, n_NaNa (r), n_NaCl (r), and n_ClCl (r) at 25^°C for the A6 MD system [c = 1.0 mol/L, each Na⁺(Cl⁻) 36 ions with 1956 water molecules] with κ = 0.25 Å⁻¹ and κ = 0.10 Å⁻¹ are shown in Figures 4 and 5, respectively. The radial distribution functions _ii (r) of other MD systems (the A3 ∼ A5 and the A7 ∼ A9) with the same κ = 0.25 Å⁻¹ are very similar to those of the A6 (Figure 4) with slightly different magnitudes of peaks, which are rather similar to the ordinary _iw (r) functions (Figure 3). However, the _ii (r) functions with different κ are totally different. For comparison, we show the _ii (r) and n_ii (r) functions of the A6 MD system with κ = 0.10 Å⁻¹ in Figure 5, in which a dramatic change occurs with a very deep minimum of _NaCl(r) and, as a result, sharp maxima of _NaNa (r) and _ClCl (r) at the distance 9.95 Å. Which is closer to nature, κ = 0.10 Å⁻¹ or 0.25 Å⁻¹?

Figure 4 

Radial distribution functions gNaNa(r), gNaCl(r), and gClCl(r), and the running coordination numbers, nNaNa(r), nNaCl(r), and nClCl(r) for the A6 MD system at 25°C as a function of r (Å) between the ions (i) using Ewald sum with κ = 0.25 Å − 1.

Figure 5 

Radial distribution functions gNaNa(r), gNaCl(r), and gClCl(r), and the running coordination numbers, nNaNa(r), nNaCl(r), and nClCl(r) for the A6 MD system at 25°C as a function of r (Å) between the ions (i) using Ewald sum with κ = 0.10 Å − 1.

Other MD systems (the A3 ∼ A5 and the A7 ∼ A9) with the same κ = 0.10 Å⁻¹ have also almost the same _ii(r) to those of the A6 (Figure 5) with slightly different magnitudes of peaks. Also, the _ii(r) functions for the A3 ∼ A9 MD systems with κ = 0.0 (no Ewald sum), 0.05 Å⁻¹, and 0.15 Å⁻¹ are similar to those for those MD system with κ = 0.10 Å⁻¹ (Figure 5), but the running coordination numbers n_ii(r) are different, while those _ii(r) for the A3 ∼ A9 MD systems with κ = 0.20 Å⁻¹ are similar to those with κ = 0.25 Å⁻¹ (Figure 4). In Figure 5, _NaNa (r) and _ClCl (r) at short distances are nonzero but the corresponding average n_NaNa (r) at 5.0 Å and n_ClCl (r) at 6.5 Å are very small, 0.17 and 0.21, respectively, but n_NaCl (r) at 5.0 Å and 6.5 Å are already 0.28 and 1.20. As the distance r increases, n_NaNa (r) and n_ClCl (r) increase slowly, then suddenly reaching 2.04 and 2.10 at 9.95 Å, while the corresponding average n_NaCl (r) is 2.73 at 8.0 Å and slowly reaches about 3.80 at 9.95 Å.

3.5. Ionic Atmosphere

We have also listed the average coordination numbers n_NaNa (R₁), n_ClCl (R₁), and n_NaCl (R₁) ( = n_ClNa (R₁)) for the A3 ∼ A9 MD systems with κ = 0.10 Å⁻¹ in the coordination shells around Na⁺ and Cl⁻ at 25°C with R₁ = 9.95 Å in Table 5, while the coordination numbers n_ii (R₁) for the A3 ∼ A9 MD systems with κ = 0.25 Å⁻¹ are not defined since _ii (r) for κ = 0.25 Å⁻¹ show neither any clear maxima for _ClCl (r) and _NaNa (r) at R₁ = 9.95 Å nor a clear minimum for _NaCl (r) at R₁ = 9.95 Å. For the A3 (each Na⁺(Cl⁻) 4 ions) and the A4 (each Na⁺(Cl⁻) 8 ions), n_NaNa (R₁) and n_ClCl (R₁) are less than 1. The ratios of n_NaCl (R₁)/(n_NaNa (R₁) + 1) and n_ClNa (R₁)/(n_ClCl (R₁) + 1) for the A5 ∼ A9 MD systems are greater than 1 which indicates a characteristic of ionic atmosphere: Averaged over time, counterions are more likely to be found by any given ion. The time averaged, spherical haze, in which counterions outnumber ions of the same charge as the central ion, has a net charge opposite in sign to that on the central ion.

Table 5 

Average coordination numbers (n_ii (R₁)) and residence times (τ_ii, ps) of Na⁺ and Cl⁻ in the first coordination shells around Na⁺ and Cl⁻ in SPC/E water at 25°C with κ = 0.10 Å⁻¹ employing a velocity Verlet algorithm where R₁ = 9.95 Å. n_ii (R1) = n_ii (R1) +1.

Relating to the running coordination numbers, n_NaNa (r), n_NaCl (r), and n_ClCl (r) as seen in Table 5, we can imagine an imaginary sphere with Na⁺ at the origin and 9.95 Å from it, filled with water molecules, with more number of Cl⁻ inside the sphere and less number of Na⁺ near the surface of the sphere. Figure 6 shows a snapshot of the ionic atmosphere with Na⁺ at the origin and 9.95 Å from it, obtained from our MD simulation. Since an electric field is not applied to a solution of ions in this study, an asymmetric ionic atmosphere, the relaxation and the electrophoretic effects are not observed, but it is currently under studying by applying an electric field to a solution of ion.

Figure 6 

A snapshot of ionic atmosphere with Na + at the origin and R1 = 9.95 Å from it, filled with 132 water molecules at 25°C. 3 Na+: dark blue and 4 Cl-: dark cyan.

3.6. Residence Time Correlation Functions and Residence times

The residence times are calculated from time correlation functions [8–10, 24] defined bywhere θ_i (r, t) is the Heaviside unit step function, which is 1 if a water molecule i is in a region r within the coordination shell around the ion at time t and 0 otherwise, and N_h is the average number of water molecules in this region r at t = 0 [8–10, 24]. Figure 7 shows the time dependence of Ln[R (r, t)] for water in the first hydration shell of the Na⁺ and Cl⁻ ions with R₁ = 3.95 Å for _ClO (r) and 3.15 Å for _NaO (r) and that for Na⁺or Cl⁻ in the coordination shell around the Cl⁻ or Na⁺ with R₁ = 9.95 Å at 25°C calculated from our MD simulations. The residence time t is obtained by fitting the time correlation function to an exponential decay <R (r, t)> exp (–t/τ), which is useful particularly when t is large. Table 4 shows the average residence times of water in the first hydration shell of the Na⁺ and Cl⁻ ions for all the MD systems with κ = 0.05 Å⁻¹, 0.10 Å⁻¹, and 0.25 Å⁻¹.

Figure 7 

Logarithm of residence time correlation functions for the hydrated SPC/E water molecules at 25°C in the first hydration shell of Na+ and Cl- and that for Na + or Cl-in the coordination shell around Cl- or Na+ with R1 = 9.95 Å obtained from our MD simulations with κ = 0.10 Å − 1.

In Table 4, the values of τ_NaO and τ_ClO for the A1 MD system are exceptionally smaller than those for the other MD systems. τ_NaO are always greater than τ_ClO with given values of κ and both the residence times increase with κ for all the systems. Roughly saying, τ_NaO with κ = 0.05 Å⁻¹ increases and decreases with the concentration of NaCl, c, and τ_ClO increases monotonically with c, while τ_NaO and τ_ClO with κ = 0.10 and 0.25 Å⁻¹ increase with c. Generally, both residence times for all the MD systems with κ = 0.05 Å⁻¹ are too low, those with κ = 0.25 Å⁻¹ are too high, and those with κ = 0.10 Å⁻¹ are moderate.

Table 5 shows the average residence times of Na⁺ and Cl⁻ in the coordination shell around Na⁺ and Cl⁻ with R₁ = 9.95 Å for the A3 ∼ A9 systems with κ = 0.10 Å⁻¹. Approximately, the values of τ_NaNa and τ_ClCl are in the range of 20–50 ps, while τ_NaCl and τ_ClNa are in that of 110 and 115 ps. Also, τ_NaNa and τ_ClCl increase with c, while τ_NaCl and τ_ClNa are independent of c. The ratios of τ_NaCl/τ_NaNa and τ_ClNa/τ_ClCl for the A5 ∼ A9 MD systems are greater than 1 which indicates that oppositely charged ions attract each other. Also, the increase of τ_NaNa and τ_ClCl and, as a result, the decrease of the ratios of τ_NaCl/τ_NaNa and τ_ClNa/τ_ClCl with c means that the migration of Na⁺(or Cl⁻) occurs hardly due to the increasing numbers of Na⁺(or Cl⁻)

3.7. Water Diffusion

We have listed diffusion coefficients D_w of SPC/E water at 25°C for the all MD systems with κ = 0.05, 0.10, and 0.25 Å⁻¹ in Table 6. For the pure water system, a recent MD simulation study [25] of 1024 SPC/E water molecules using Ewald with κ = 0.20 Å⁻¹ employing the vV algorithm reported 2.78 × 10⁻⁹ m²/s obtained from MSD (mean square displacement, equation (9)) and 2.75 × 10⁻⁹ m²/s from VAC (velocity autocorrelation, equation (10)) function at 300 K, which overestimate the experimental data at the same temperature (2.39 × 10⁻⁹ m²/s [26, 27] and 2.49 × 10⁻⁹ m²/s [26–28]). However, TIP4P/2005 water model [29] gives an excellent agreement, 2.39 × 10⁻⁹ m²/s, obtained from both the MSD and VAC in our previous study [25] with the experiment one.

Table 6 

Water diffusion coefficients (10⁻⁹ m²/s) at 25°C employing a velocity Verlet algorithm, where n is the number of (Na⁺ + Cl⁻) ions.

D_w obtained from our recent MD simulations [8] as a function of RFC(reaction field correction), ST(simple truncation), and Ewald sum parameter employing the Gear algorithm overestimates the experimental data, while D_w employing the vV algorithm underestimates except for ST and Ewald with κ > 0.45 Å⁻¹. These MD simulations [8] reported that D_w is 3.03 × 10⁻⁹ m²/s in the Na⁺-water system and 3.09 × 10⁻⁹ m²/s in the Cl⁻-water system using ST employing the Gear algorithm with N_w = 215 of SPC/E water. An MD simulation for alkali Earth metal cations (Mg²⁺, Ca²⁺, Sr²⁺, and Ba²⁺) in an aqueous solution at 25^°C of SPC/E water potential (N_w = 215) using Ewald employing the Gear algorithm reported that D_w is in the range of 2.46–2.62 × 10⁻⁹ m²/s [17].

In Table 6, roughly speaking, D_w with κ = 0.05 Å⁻¹ increases up to the A7 MD system with slow decreases as c increases, while D_w with κ = 0.10 Å⁻¹ has almost the same values up to the A6 MD system with sudden decreases, and D_w with κ = 0.25 Å⁻¹ decreases monotonically up to the A6 MD system with sudden decreases. The behavior of D_w with κ = 0.10 Å⁻¹ seems to be the most reasonable one: D_w has almost the same values up to n = 36 Na⁺(Cl⁻) ions and decreases suddenly for n = 72, 108, and 144.

In a study of the system-size dependence of translational diffusion coefficients and viscosities of pure water in molecular dynamics simulations under periodic boundary conditions [30], for a cubic simulation box of length L, the diffusion coefficient is corrected for system-size effects as D₀ = D_PBC + 2.837297 k_BT/(6πηL), where D_PBC is the diffusion coefficient calculated in the simulation, k_B is the Boltzmann constant, T is the absolute temperature, and η is the shear viscosity of the solvent. Using T = 298.15 K and L = 39.1009 Å for our systems and η = 0.854 cP at 300 K, D₀-D_PBC = 0.1856 × 10⁻⁹ m²/s which is quite small, compared to our calculated D_w = 2.02–2.84 × 10⁻⁹ m²/s with κ = 0.10 Å⁻¹ (Table 6). This correction for D_w is not included in this study.

A recent study [31] reported that the effect of salt on the dynamics of water molecules follows the Hofmeister series. For some structure-making salts, the self-diffusion coefficient of the water molecules, D_w, decreases with increasing salt concentration and for other structure-breaking salts, D_w increases with increasing salt concentration c. Both ratios of experimental and simulation D_w (NaCl)/D₀ (bulk) almost linearly decrease with c up to 0.75 and 0.5 at c = 4 M.

A more recent study [32] reported that when used with a good-quality water model, e.g., TIP4P/2005 [29] or E3B [33], this method recovers the qualitative behavior of the water diffusion trends of experiment shown that the water diffusion coefficient increases in the presence of salts of low charge density (e.g., CsI), whereas the results of simulations with nonpolarizable models show a decrease of the water diffusion coefficient in all alkali halide solutions.

3.8. Energetics

Equations (7) and (8) have 18 potential energies, which are (1) water-water, (2) water-Na⁺, (3) water-Cl^-, (4) Na⁺-Na⁺, (5) Cl⁻-Cl^-, and (6) Na⁺-Cl^- LJ potential energy and the corresponding electrostatic (Coulomb) energy in real and reciprocal spaces for the Ewald sum. The important potential energy that determines the molar conductivity of an electrolyte can be the electrostatic energy of ion-ion in real space and reciprocal spaces for the Ewald sum. In Table 7, we compared the average electrostatic energies of the electrolytic system obtained using the Ewald sum of κ = 0.10 Å⁻¹ and the velocity Verlet algorithm. Na⁺-Na⁺ and Cl⁻-Cl^- electrostatic energies in real space increase monotonically as the number of Na⁺(Cl⁻) ions increases and the corresponding Na⁺-Cl^- energy also increases negatively, whereas the corresponding electrostatic energies in reciprocal space increase and decrease with the number of Na⁺(Cl⁻) ions. The ion-ion electrostatic energy in real space seems to be more important than in reciprocal space when determining the molar conductivity of an electrolyte. Ion-ion electrostatic energy, which increases monotonically in real space with the number of Na⁺(Cl⁻) ions, positively or negatively, retards the diffusion of ions. In conclusion, the molar conductivity of the electrolyte decreases with the concentration of the electrolyte.

Table 7 

The average electrostatic energies (kJ/mol) per ion of the electrolyte systems using Ewald sum with κ = 0.10 Å⁻¹ employing a velocity Verlet algorithm, where n is the number of (Na⁺ or Cl⁻) ion(s).

The molar conductivity of an electrolyte is found to vary with the concentration. One reason for this variation is that the number of ions in the solution might not be proportional to the concentration of the electrolyte (that is, a strong or a weak electrolyte). Second, because ions interact with each other strongly, the conductivity of a solution is not exactly proportional to the number of ions present. Figure 8 shows the sum of diffusion coefficients (○) of all ions multiplied by n in D_A1 to D_A9 in Table 3 for systems A1 to A9 for κ = 0.10 Å⁻¹. The sum of the “Ideal” diffusion coefficients (---) of all ions obtained from D_A1 multiplied by n in Table 3 for all identical systems is also shown. Here, of course, the molar conductivity of a solution is exactly proportional to the number of ions present. The deviation of the sum (○) of the diffusion coefficients from the “Ideal” diffusion coefficient (---) is due to the delay of the moving ions due to the strong ion-ion electrostatic interaction in real space as described above.

Figure 8 

Sum of diffusion coefficients of all the ions in the systems. MD result: multiplied by n in DA1 to DA9 in Table 3 for systems A1 to A9 for κ = 0.10 Å − 1 and ideal: obtained from DA1 multiplied by n in Table 3 for all identical systems.

4. Conclusions

We have carried out molecular dynamics (MD) simulations of NaCl in SPC/E water at 25^°C to calculate the molar conductivity of NaCl as a function of NaCl concentration (c) using Ewald sums employing a velocity Verlet algorithm.

It is found that the structural properties of the MD systems with different κ are totally different. The obtained ion-ion radial distribution functions _ii (r) with κ = 0.25 Å⁻¹ show the ordinary behavior similar to the ion-water radial distribution functions _iw (r). However, those functions _ii(r) with κ = 0.10 Å⁻¹ give a dramatic change with a very deep minimum of _NaCl(r) and, as a result, sharp maxima of _NaNa (r) and _ClCl (r) at the distance R₁ = 9.95 Å. The ratios of n_NaCl (R₁)/(n_NaNa (R₁) + 1) and n_ClNa (R₁)/(n_ClCl (R₁)+1) for the A5 ∼ A9 MD systems with κ = 0.10 Å⁻¹ are greater than 1 which indicates a characteristic of ionic atmosphere. This is the actual evidence of the basis of the Debye–Hückel theory of ionic solutions.

It is also found that the MD result for Λ_m of NaCl with Ewald sum parameter κ = 0.10 Å⁻¹ gives the closest one to the experimental data. The behavior of is close to the experimental date except at very low concentrations which are much higher than the experimental data. The analysis of radial distribution functions, hydration numbers, coordination numbers around Na⁺ and Cl⁻, residence times of water around Na⁺ and Cl⁻, water diffusion, and ion-ion electrostatic energies support this finding. It has been found that the ion-ion electrostatic energy in real space plays an important role in explaining that the molar conductivity of the electrolyte decreases with concentration.

Compared with the experimental data, the best MD results of and for N_NaCl = 1, 2, 4, 8, and 18 are obtained with κ = 0.25 Å⁻¹. This is consistent with the best result for µ_Na+ and µ_Cl- in single-ion systems in 1023 SPC/E water molecules. [2] This indicates that the Ewald sum parameter κ = 0.25 Å⁻¹ is suitable for systems with a small number of ions and 1000 to 2000 water molecules. However, as the number of ions increases without significant fluctuations in the number of water molecules, κ changes for the best MD results of and . We conclude that the appropriate value of κ is determined by the number of ions rather than the water molecule.

Data Availability

The data used in this study are given in the manuscript.

Conflicts of Interest

The author declares no conflicts of interest.