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

Authors: Song Hi Lee
Journal: Journal of Chemistry
 (2020)

---

## Abstract

We present the result of molecular dynamics (MD) simulations to calculate the molar conductivityΛm=λNa++λCl− 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 Å−1 gives the closest one to the experimental data and that the obtained radial distribution functions gii(r) with κ = 0.10 Å−1 show a dramatic change with a very deep minimum of gNaCl(r) and, as a result, sharp maxima of gNaNa(r) and gClCl(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.


---

## Body

## 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 Λm 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 Λmo in an infinitely dilute solution is related to the limiting molar conductivity (λ+o and λ−o) of each ion, as shown in equation (2) (Kohlrausch’s law of the independent migration of ions).(1)Λm=Λmo−Κc1/2,where K is a constant.(2)Λmo=ν+λ+o+ν−λ−o,where the limiting molar conductivities of Na+ and Cl− at 25oC are measured as 5.01 and 7.63 mS m2 mol−1 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:(3)Λm=Λmo−A+BΛmoc1/2,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 m2 mol−1/(mol·L−1)1/2 and 0.229 (mol L−1)1/2 [1]A systematic conductivity study of sodium chloride-water + 1, 4-dioxane in dilute solutions (c < 0.01 mol-dm−3) [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 KA 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−3 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 Å−1 employing a velocity Verlet algorithm gives the best result for μNa+ and μCl− at 25 C, even though the calculated μNa+ overestimates the experimental data and the μCl− underestimates them. This behavior of μNa+ 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 ofc = 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 = Vsoln = 36000 cm3/NA with NA, 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/cm3) of NaCl solution in the literature. Instead, we use the experimental data of the weight percent (%WNaCl) and the density (ρt, g/cm3) of NaCl solutions (Table 1) to find the numbers of ions (NNaCl) and water (Nw). First, the relationship between ρt and %WNaCl is obtained by using the least squares method from Table 1.(4)ρt=0.00758333853%WNaCl+0.9942949g/cm3.Table 1
Experimental data between %WNaCl and ρt at 25oC. The relationship (ρt = 0.00758333853 %WNaCl + 0.9942949) is obtained by using the least squares method.
%WNaCl124681012ρt (g/cm3)1.004091.011121.025301.039631.054121.068791.08356%WNaCl14161820222426ρt (g/cm3)1.098721.114011.129541.145331.161401.177601.19443Forc = 1 mol/L of NNaCl = 36 with V = 36000 cm3/NA,(5)ρt=18.0152Nw+36×58.443/36000g/cm3.Also,(6)%WNaCl=36×58.443/18.0152Nw+36×58.443×100.Substituting equations (5) and (6) into equation (4), we obtain Nw = 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.
SystemNumber of Na+(Cl−) ion (s)Number of molecules (Nw)Number of water (mol/L)%WNaClρt (g/cm3)A1119860.02780.163080.99546A2219850.05560.325800.99658A3419840.1110.649800.99933A4819810.2221.293141.00432A51819720.52.875981.01605A63619561.05.634321.03727A77219182.010.85601.07670A810818753.015.74411.11362A914418284.020.35381.14854The SPC/E model [10] was adopted for water-water and ion-water. The pair potential between water and ion has the TIPS form [11]:(7)viw=4εioσiorio12−σiorio6+∑jεwqiqjrij,where σio and εio are Lennard–Jones (LJ) parameters between oxygen on a water molecule and an ion i, qj is the charge at site j in water, and qi is the charge on ion i. Moreover, rio and rij 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 qO2−=−0.8476e, qH+=0.4238e, qNa+=0.1e, and qCl−=−0.1e. The pair potential energy of ion-ion also has the TIPS form [11]:(8)vij=4εijσijrij12−σijrij6+12∑j≠iqiqjrij,where σij and εij are Lennard–Jones (LJ) parameters between ion i and ion j, qi is the charge on ion i, and rij 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 Table1, 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 Å−1 in 0.05 Å−1 increments, and the cutoff sphere of radius Rc in the real space and kmax in the reciprocal space are chosen as L/2 and 7, but |k|2 ≤ 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),(9)Di=16limt⟶∞ddt<rit−ri02>,and the Green–Kubo (GK) relation from the velocity autocorrelation (VAC) function,(10)Di=13∫0∞dt<vi0⋅vit>.


## 3. Results and Discussion

### 3.1. Diffusion Coefficient

DNa+ and DCl− 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 Å−1 employing a velocity Verlet algorithm. We have listed sum of Na+ and Cl− diffusion coefficients D=DNa++DCl− 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 Å−1 at 3.0 and 4.0 mol/L and those (□) with κ = 0.05 Å−1 at 2.0, 3.0, and 4.0 mol/L.Table 3
Diffusion coefficients (D, 10−9 m2·sec−1) 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).
System (n)κ = 0.05κ = 0.10κ = 0.15κ = 0.20κ = 0.25Best to exp.(κ, Å−1)A1 (1)5.735.065.084.564.61A1 (0.25)4.61A2 (2)4.653.853.823.943.73A2 (0.25)3.73A3 (4)3.913.463.643.473.15A3 (0.25)3.15A4 (8)3.503.203.283.222.98A4 (0.25)2.98A5 (18)3.052.823.173.002.69A5 (0.25)2.69A6 (36)2.432.482.862.862.65A6 (0.05)2.43A7 (72)1.812.052.462.482.36A7 (0.10)2.05A8 (108)1.411.642.072.102.02A8 (0.10)1.64A9 (144)1.141.341.691.781.73A9 (0.10)1.34Figure 1
Comparison of the diffusion coefficients (D =  DNa+ + DCl−) 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 forD are obtained at 0.0278–0.50 mol/L with κ = 0.25 Å−1 (▽), at 1.0 M with κ = 0.05 Å−1 (□), and at 2.0–4.0 mol/L with κ = 0.10 Å−1 (◇). 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 Å−1 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 Nw = 1800–2000. Note that the best results of D for NNaCl = 1, 2, 4, 8, and 18 (the A1 ∼ A5 MD systems) are obtained with κ = 0.25 Å−1 which also gives the best result for µNa+ and µCl- in single-ion systems in 1023 SPC/E water molecules. [2].In Figure2, we have compared the best obtained D (□) and D with κ = 0.10 Å−1 (◇) 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 Å−1 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 DNa+ and DCl− with κ = 0.10 Å−1 as a function of c. DCl- at all the c (mol/L) is greater than DNa+, and these two values become equal as c increases.Figure 2
Comparison of the best diffusion coefficients (D), those withκ = 0.10 Å-1, and DNa+ and DCl− 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, theD 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Λm=λNa++λCl− of NaCl as a function of c1/2 obtained from D in our MD simulations in the inset of Figure 2, where λi=Dizi2F2/RT, where zi 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Λm of NaCl show two linear dependences (Ex1 and Ex2) as a function of c1/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 m2 mol−1/(mol·L−1)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 m2 mol−1/(mol·L−1)1/2, respectively. According to the Debye–Huckel–Onsager (DHO) theory, equation (3), the slope A+BΛmo is equal to -8.91 mS m2 mol−1/(mol·L−1)1/2 with the numerical values of A and B in water, 6.02 mS m2 mol−1/(mol·L−1)1/2 and 0.229 (mol L−1)1/2, and Λmo of NaCl, 12.64 mS m2 mol−1 [1]. In conclusion, the MD results for Λm 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 Functionsgiwr and Running Hydration Numbers niw (r)

The radial distribution functions,gio (r) and gih (r), and the running hydration numbers, nio(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 Å−1 as a function of distance. The running hydration number is defined as(11)nijr=∫0rgijr′4πr'2dr′.Figure 3
Radial distribution functionsgio(r) and gih(r) and the running hydration numbers nio (r) and nih (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.Figure3 shows the typical giw (r) functions with positions of gClO (r) at the first maximum 3.25 Å and the first minimum 3.95 Å, of gCl-H (r) at the first maximum 2.30 Å and the first minimum 3.10 Å, gNaO (r) at the first maximum 2.35 Å and the first minimum 3.15 Å and gNaH (r) at the first maximum 3.00 Å and the first minimum 3.75 Å. The positions of the first maxima and minima of gio (r) and gih (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 giw (r) of multi-ion systems are almost the same to those of single-ion systems [2].The hydration numbernio (R1) in the first shell is defined through equation.(11) from the ion-oxygen distribution functions gio (r) where the upper limit of integration R1 is the radius of the first hydration sphere, which corresponds to the first minimum in gio (r). That is, R1 = 3.95 Å for gClO (r) and 3.15 Å for gNaO (r). We have listed the average hydration numbers nNaO (R1) and nClO (R1) for the values of κ = 0.05, 0.10 and 0.25 Å−1 in the first hydration shells of Na+ and Cl− in SPC/E water at 25°C in Table 4. The values of nNaO (R1) and nClO (R1) are almost the same for κ = 0.05 Å−1 (5.64–5.85 and 7.55–7.62), κ = 0.10 Å−1 (5.67–5.90 and 7.65–7.82), and κ = 0.25 Å−1 (5.47–5.69 and 7.51–7.68).Table 4
The average hydration numbers (nio (R1)) 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 R1 = 3.15 Å for Na+ and R1 = 3.95 Å for Cl−.
Systemκ = 0.05 Å−1κ = 0.10 Å−1κ = 0.25 Å−1nNaO (R1)τNaOnClO (R1)τClOnNaO (R1)τNaOnClO (R1)τClOnNaO (R1)τNaOnClO (R1)τClOA15.6512.97.628.595.6713.57.829.265.6914.07.589.26A25.6615.67.639.525.7018.07.829.865.6518.67.6410.4A35.6616.77.639.385.7118.77.8210.35.6819.47.5112.0A45.6616.37.599.565.7018.27.8110.65.6921.37.5212.1A55.6414.37.579.435.7018.47.7810.55.6721.27.5412.8A65.6713.87.559.475.7418.07.7811.15.6422.47.5913.1A75.7813.27.599.755.7617.97.6511.75.5323.77.6314.2A85.8514.07.6210.35.9018.87.6912.65.4725.87.6815.6A95.9413.97.6510.85.9820.07.7214.15.3128.07.7117.1


### 3.4. Radial Distribution Functionsgii (r) and Running Coordination Numbers nii (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 gii(r) functions due to the very small numbers of ion(s). The radial distribution functions, gNaNa (r), gNaCl (r), and gClCl (r), and the running coordination numbers, nNaNa (r), nNaCl (r), and nClCl (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 Å−1 and κ = 0.10 Å−1 are shown in Figures 4 and 5, respectively. The radial distribution functions gii (r) of other MD systems (the A3 ∼ A5 and the A7 ∼ A9) with the same κ = 0.25 Å−1 are very similar to those of the A6 (Figure 4) with slightly different magnitudes of peaks, which are rather similar to the ordinary giw (r) functions (Figure 3). However, the gii (r) functions with different κ are totally different. For comparison, we show the gii (r) and nii (r) functions of the A6 MD system with κ = 0.10 Å−1 in Figure 5, in which a dramatic change occurs with a very deep minimum of gNaCl(r) and, as a result, sharp maxima of gNaNa (r) and gClCl (r) at the distance 9.95 Å. Which is closer to nature, κ = 0.10 Å−1 or 0.25 Å−1?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 Å−1 have also almost the same gii(r) to those of the A6 (Figure 5) with slightly different magnitudes of peaks. Also, the gii(r) functions for the A3 ∼ A9 MD systems with κ = 0.0 (no Ewald sum), 0.05 Å−1, and 0.15 Å−1 are similar to those for those MD system with κ = 0.10 Å−1 (Figure 5), but the running coordination numbers nii(r) are different, while those gii(r) for the A3 ∼ A9 MD systems with κ = 0.20 Å−1 are similar to those with κ = 0.25 Å−1 (Figure 4). In Figure 5, gNaNa (r) and gClCl (r) at short distances are nonzero but the corresponding average nNaNa (r) at 5.0 Å and nClCl (r) at 6.5 Å are very small, 0.17 and 0.21, respectively, but nNaCl (r) at 5.0 Å and 6.5 Å are already 0.28 and 1.20. As the distance r increases, nNaNa (r) and nClCl (r) increase slowly, then suddenly reaching 2.04 and 2.10 at 9.95 Å, while the corresponding average nNaCl (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 numbersnNaNa (R1), nClCl (R1), and nNaCl (R1) ( = nClNa (R1)) for the A3 ∼ A9 MD systems with κ = 0.10 Å−1 in the coordination shells around Na+ and Cl− at 25°C with R1 = 9.95 Å in Table 5, while the coordination numbers nii (R1) for the A3 ∼ A9 MD systems with κ = 0.25 Å−1 are not defined since gii (r) for κ = 0.25 Å−1 show neither any clear maxima for gClCl (r) and gNaNa (r) at R1 = 9.95 Å nor a clear minimum for gNaCl (r) at R1 = 9.95 Å. For the A3 (each Na+(Cl−) 4 ions) and the A4 (each Na+(Cl−) 8 ions), nNaNa (R1) and nClCl (R1) are less than 1. The ratios of nNaCl (R1)/(nNaNa (R1) + 1) and nClNa (R1)/(nClCl (R1) + 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 (nii (R1)) 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 Å−1 employing a velocity Verlet algorithm where R1 = 9.95 Å. nii (R1)∗ = nii (R1) +1.
SystemnNaNa (R1)τNaNanClCl (R1)τClClnNaC l (R1)τNaClnClNa (R1)τClNanNaCl (R1)/nNaNa (R1)∗nClNa (R1)/nClCl (R1)∗τNaCl/τNaNaτClNa/τClClA3.11812.8.18214.6.63255.3.63258.0——4.323.97A4.37115.2.37417.61.2274.51.2294.0——4.905.34A51.2522.11.4125.02.70110.2.70120.1.201.124.984.80A62.0427.02.1029.63.80115.3.80114.1.251.234.263.85A73.6734.33.9538.15.84102.5.84102.1.251.182.972.68A85.8342.26.3348.28.29106.8.29108.1.211.132.512.24A97.9348.68.5551.310.7115.10.7115.1.201.122.372.24Relating to the running coordination numbers,nNaNa (r), nNaCl (r), and nClCl (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 andR1 = 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 by(12)Rr,t=1Nh∑i=1Nhθir,0⋅θir,t,where θ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 Nh 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 R1 = 3.95 Å for gClO (r) and 3.15 Å for gNaO (r) and that for Na+or Cl− in the coordination shell around the Cl− or Na+ with R1 = 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 Å−1, 0.10 Å−1, and 0.25 Å−1.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+ withR1 = 9.95 Å obtained from our MD simulations with κ = 0.10 Å − 1.In Table4, 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 Å−1 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 Å−1 increase with c. Generally, both residence times for all the MD systems with κ = 0.05 Å−1 are too low, those with κ = 0.25 Å−1 are too high, and those with κ = 0.10 Å−1 are moderate.Table5 shows the average residence times of Na+ and Cl− in the coordination shell around Na+ and Cl− with R1 = 9.95 Å for the A3 ∼ A9 systems with κ = 0.10 Å−1. 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 coefficientsDw of SPC/E water at 25°C for the all MD systems with κ = 0.05, 0.10, and 0.25 Å−1 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 Å−1 employing the vV algorithm reported 2.78 × 10−9 m2/s obtained from MSD (mean square displacement, equation (9)) and 2.75 × 10−9 m2/s from VAC (velocity autocorrelation, equation (10)) function at 300 K, which overestimate the experimental data at the same temperature (2.39 × 10−9 m2/s [26, 27] and 2.49 × 10−9 m2/s [26–28]). However, TIP4P/2005 water model [29] gives an excellent agreement, 2.39 × 10−9 m2/s, obtained from both the MSD and VAC in our previous study [25] with the experiment one.Table 6
Water diffusion coefficients (10−9 m2/s) at 25°C employing a velocity Verlet algorithm, where n is the number of (Na+ + Cl−) ions.
System (n)κ = 0.05 Å−1κ = 0.10 Å−1κ = 0.25 Å−1A1 (2)2.742.812.68A2 (4)2.762.792.66A3 (8)2.802.822.64A4 (16)2.872.822.62A5 (36)3.012.842.52A6 (72)3.152.802.38A7 (144)3.252.582.05A8 (216)3.142.321.76A9 (288)2.942.021.51Dw 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 Dw employing the vV algorithm underestimates except for ST and Ewald with κ > 0.45 Å−1. These MD simulations [8] reported that Dw is 3.03 × 10−9 m2/s in the Na+-water system and 3.09 × 10−9 m2/s in the Cl−-water system using ST employing the Gear algorithm with Nw = 215 of SPC/E water. An MD simulation for alkali Earth metal cations (Mg2+, Ca2+, Sr2+, and Ba2+) in an aqueous solution at 25°C of SPC/E water potential (Nw = 215) using Ewald employing the Gear algorithm reported that Dw is in the range of 2.46–2.62 × 10−9 m2/s [17].In Table6, roughly speaking, Dw with κ = 0.05 Å−1 increases up to the A7 MD system with slow decreases as c increases, while Dw with κ = 0.10 Å−1 has almost the same values up to the A6 MD system with sudden decreases, and Dw with κ = 0.25 Å−1 decreases monotonically up to the A6 MD system with sudden decreases. The behavior of Dw with κ = 0.10 Å−1 seems to be the most reasonable one: Dw 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 D0 = DPBC + 2.837297 kBT/(6πηL), where DPBC is the diffusion coefficient calculated in the simulation, kB 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, D0-DPBC = 0.1856 × 10−9 m2/s which is quite small, compared to our calculated Dw = 2.02–2.84 × 10−9 m2/s with κ = 0.10 Å−1 (Table 6). This correction for Dw 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, Dw, decreases with increasing salt concentration and for other structure-breaking salts, Dw increases with increasing salt concentration c. Both ratios of experimental and simulation Dw (NaCl)/D0 (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 Å−1 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 Å−1 employing a velocity Verlet algorithm, where n is the number of (Na+ or Cl−) ion(s).
System (n)Na+-Na+Cl−-Cl-Na+-Cl-Real reciprocalReal reciprocalReal reciprocalA1 (1)————−0.293.83A2 (2)0.0013.40.9315.7−11.0−2.16A3 (4)1.7030.92.8926.5−23.3−7.59A4 (8)6.8633.26.6734.9−48.2−18.0A5 (18)34.138.029.939.3−100−42.8A6 (36)52.834.542.731.7−131−32.9A7 (72)10715.990.214.1−193−13.5A8 (108)18214.416212.8−272−11.9A9 (144)2574.912324.61−347−4.05The 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. Figure8 shows the sum of diffusion coefficients (○) of all ions multiplied by n in DA1 to DA9 in Table 3 for systems A1 to A9 for κ = 0.10 Å−1. The sum of the “Ideal” diffusion coefficients (---) of all ions obtained from DA1 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 Table3 for systems A1 to A9 for κ = 0.10 Å − 1 and ideal: obtained from DA1 multiplied by n in Table 3 for all identical systems.


## 3.1. Diffusion Coefficient

DNa+ and DCl− 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 Å−1 employing a velocity Verlet algorithm. We have listed sum of Na+ and Cl− diffusion coefficients D=DNa++DCl− 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 Å−1 at 3.0 and 4.0 mol/L and those (□) with κ = 0.05 Å−1 at 2.0, 3.0, and 4.0 mol/L.Table 3
Diffusion coefficients (D, 10−9 m2·sec−1) 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).
System (n)κ = 0.05κ = 0.10κ = 0.15κ = 0.20κ = 0.25Best to exp.(κ, Å−1)A1 (1)5.735.065.084.564.61A1 (0.25)4.61A2 (2)4.653.853.823.943.73A2 (0.25)3.73A3 (4)3.913.463.643.473.15A3 (0.25)3.15A4 (8)3.503.203.283.222.98A4 (0.25)2.98A5 (18)3.052.823.173.002.69A5 (0.25)2.69A6 (36)2.432.482.862.862.65A6 (0.05)2.43A7 (72)1.812.052.462.482.36A7 (0.10)2.05A8 (108)1.411.642.072.102.02A8 (0.10)1.64A9 (144)1.141.341.691.781.73A9 (0.10)1.34Figure 1
Comparison of the diffusion coefficients (D =  DNa+ + DCl−) 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 forD are obtained at 0.0278–0.50 mol/L with κ = 0.25 Å−1 (▽), at 1.0 M with κ = 0.05 Å−1 (□), and at 2.0–4.0 mol/L with κ = 0.10 Å−1 (◇). 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 Å−1 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 Nw = 1800–2000. Note that the best results of D for NNaCl = 1, 2, 4, 8, and 18 (the A1 ∼ A5 MD systems) are obtained with κ = 0.25 Å−1 which also gives the best result for µNa+ and µCl- in single-ion systems in 1023 SPC/E water molecules. [2].In Figure2, we have compared the best obtained D (□) and D with κ = 0.10 Å−1 (◇) 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 Å−1 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 DNa+ and DCl− with κ = 0.10 Å−1 as a function of c. DCl- at all the c (mol/L) is greater than DNa+, and these two values become equal as c increases.Figure 2
Comparison of the best diffusion coefficients (D), those withκ = 0.10 Å-1, and DNa+ and DCl− 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, theD 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Λm=λNa++λCl− of NaCl as a function of c1/2 obtained from D in our MD simulations in the inset of Figure 2, where λi=Dizi2F2/RT, where zi 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Λm of NaCl show two linear dependences (Ex1 and Ex2) as a function of c1/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 m2 mol−1/(mol·L−1)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 m2 mol−1/(mol·L−1)1/2, respectively. According to the Debye–Huckel–Onsager (DHO) theory, equation (3), the slope A+BΛmo is equal to -8.91 mS m2 mol−1/(mol·L−1)1/2 with the numerical values of A and B in water, 6.02 mS m2 mol−1/(mol·L−1)1/2 and 0.229 (mol L−1)1/2, and Λmo of NaCl, 12.64 mS m2 mol−1 [1]. In conclusion, the MD results for Λm 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 Functionsgiwr and Running Hydration Numbers niw (r)

The radial distribution functions,gio (r) and gih (r), and the running hydration numbers, nio(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 Å−1 as a function of distance. The running hydration number is defined as(11)nijr=∫0rgijr′4πr'2dr′.Figure 3
Radial distribution functionsgio(r) and gih(r) and the running hydration numbers nio (r) and nih (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.Figure3 shows the typical giw (r) functions with positions of gClO (r) at the first maximum 3.25 Å and the first minimum 3.95 Å, of gCl-H (r) at the first maximum 2.30 Å and the first minimum 3.10 Å, gNaO (r) at the first maximum 2.35 Å and the first minimum 3.15 Å and gNaH (r) at the first maximum 3.00 Å and the first minimum 3.75 Å. The positions of the first maxima and minima of gio (r) and gih (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 giw (r) of multi-ion systems are almost the same to those of single-ion systems [2].The hydration numbernio (R1) in the first shell is defined through equation.(11) from the ion-oxygen distribution functions gio (r) where the upper limit of integration R1 is the radius of the first hydration sphere, which corresponds to the first minimum in gio (r). That is, R1 = 3.95 Å for gClO (r) and 3.15 Å for gNaO (r). We have listed the average hydration numbers nNaO (R1) and nClO (R1) for the values of κ = 0.05, 0.10 and 0.25 Å−1 in the first hydration shells of Na+ and Cl− in SPC/E water at 25°C in Table 4. The values of nNaO (R1) and nClO (R1) are almost the same for κ = 0.05 Å−1 (5.64–5.85 and 7.55–7.62), κ = 0.10 Å−1 (5.67–5.90 and 7.65–7.82), and κ = 0.25 Å−1 (5.47–5.69 and 7.51–7.68).Table 4
The average hydration numbers (nio (R1)) 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 R1 = 3.15 Å for Na+ and R1 = 3.95 Å for Cl−.
Systemκ = 0.05 Å−1κ = 0.10 Å−1κ = 0.25 Å−1nNaO (R1)τNaOnClO (R1)τClOnNaO (R1)τNaOnClO (R1)τClOnNaO (R1)τNaOnClO (R1)τClOA15.6512.97.628.595.6713.57.829.265.6914.07.589.26A25.6615.67.639.525.7018.07.829.865.6518.67.6410.4A35.6616.77.639.385.7118.77.8210.35.6819.47.5112.0A45.6616.37.599.565.7018.27.8110.65.6921.37.5212.1A55.6414.37.579.435.7018.47.7810.55.6721.27.5412.8A65.6713.87.559.475.7418.07.7811.15.6422.47.5913.1A75.7813.27.599.755.7617.97.6511.75.5323.77.6314.2A85.8514.07.6210.35.9018.87.6912.65.4725.87.6815.6A95.9413.97.6510.85.9820.07.7214.15.3128.07.7117.1


## 3.4. Radial Distribution Functionsgii (r) and Running Coordination Numbers nii (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 gii(r) functions due to the very small numbers of ion(s). The radial distribution functions, gNaNa (r), gNaCl (r), and gClCl (r), and the running coordination numbers, nNaNa (r), nNaCl (r), and nClCl (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 Å−1 and κ = 0.10 Å−1 are shown in Figures 4 and 5, respectively. The radial distribution functions gii (r) of other MD systems (the A3 ∼ A5 and the A7 ∼ A9) with the same κ = 0.25 Å−1 are very similar to those of the A6 (Figure 4) with slightly different magnitudes of peaks, which are rather similar to the ordinary giw (r) functions (Figure 3). However, the gii (r) functions with different κ are totally different. For comparison, we show the gii (r) and nii (r) functions of the A6 MD system with κ = 0.10 Å−1 in Figure 5, in which a dramatic change occurs with a very deep minimum of gNaCl(r) and, as a result, sharp maxima of gNaNa (r) and gClCl (r) at the distance 9.95 Å. Which is closer to nature, κ = 0.10 Å−1 or 0.25 Å−1?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 Å−1 have also almost the same gii(r) to those of the A6 (Figure 5) with slightly different magnitudes of peaks. Also, the gii(r) functions for the A3 ∼ A9 MD systems with κ = 0.0 (no Ewald sum), 0.05 Å−1, and 0.15 Å−1 are similar to those for those MD system with κ = 0.10 Å−1 (Figure 5), but the running coordination numbers nii(r) are different, while those gii(r) for the A3 ∼ A9 MD systems with κ = 0.20 Å−1 are similar to those with κ = 0.25 Å−1 (Figure 4). In Figure 5, gNaNa (r) and gClCl (r) at short distances are nonzero but the corresponding average nNaNa (r) at 5.0 Å and nClCl (r) at 6.5 Å are very small, 0.17 and 0.21, respectively, but nNaCl (r) at 5.0 Å and 6.5 Å are already 0.28 and 1.20. As the distance r increases, nNaNa (r) and nClCl (r) increase slowly, then suddenly reaching 2.04 and 2.10 at 9.95 Å, while the corresponding average nNaCl (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 numbersnNaNa (R1), nClCl (R1), and nNaCl (R1) ( = nClNa (R1)) for the A3 ∼ A9 MD systems with κ = 0.10 Å−1 in the coordination shells around Na+ and Cl− at 25°C with R1 = 9.95 Å in Table 5, while the coordination numbers nii (R1) for the A3 ∼ A9 MD systems with κ = 0.25 Å−1 are not defined since gii (r) for κ = 0.25 Å−1 show neither any clear maxima for gClCl (r) and gNaNa (r) at R1 = 9.95 Å nor a clear minimum for gNaCl (r) at R1 = 9.95 Å. For the A3 (each Na+(Cl−) 4 ions) and the A4 (each Na+(Cl−) 8 ions), nNaNa (R1) and nClCl (R1) are less than 1. The ratios of nNaCl (R1)/(nNaNa (R1) + 1) and nClNa (R1)/(nClCl (R1) + 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 (nii (R1)) 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 Å−1 employing a velocity Verlet algorithm where R1 = 9.95 Å. nii (R1)∗ = nii (R1) +1.
SystemnNaNa (R1)τNaNanClCl (R1)τClClnNaC l (R1)τNaClnClNa (R1)τClNanNaCl (R1)/nNaNa (R1)∗nClNa (R1)/nClCl (R1)∗τNaCl/τNaNaτClNa/τClClA3.11812.8.18214.6.63255.3.63258.0——4.323.97A4.37115.2.37417.61.2274.51.2294.0——4.905.34A51.2522.11.4125.02.70110.2.70120.1.201.124.984.80A62.0427.02.1029.63.80115.3.80114.1.251.234.263.85A73.6734.33.9538.15.84102.5.84102.1.251.182.972.68A85.8342.26.3348.28.29106.8.29108.1.211.132.512.24A97.9348.68.5551.310.7115.10.7115.1.201.122.372.24Relating to the running coordination numbers,nNaNa (r), nNaCl (r), and nClCl (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 andR1 = 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 by(12)Rr,t=1Nh∑i=1Nhθir,0⋅θir,t,where θ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 Nh 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 R1 = 3.95 Å for gClO (r) and 3.15 Å for gNaO (r) and that for Na+or Cl− in the coordination shell around the Cl− or Na+ with R1 = 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 Å−1, 0.10 Å−1, and 0.25 Å−1.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+ withR1 = 9.95 Å obtained from our MD simulations with κ = 0.10 Å − 1.In Table4, 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 Å−1 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 Å−1 increase with c. Generally, both residence times for all the MD systems with κ = 0.05 Å−1 are too low, those with κ = 0.25 Å−1 are too high, and those with κ = 0.10 Å−1 are moderate.Table5 shows the average residence times of Na+ and Cl− in the coordination shell around Na+ and Cl− with R1 = 9.95 Å for the A3 ∼ A9 systems with κ = 0.10 Å−1. 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 coefficientsDw of SPC/E water at 25°C for the all MD systems with κ = 0.05, 0.10, and 0.25 Å−1 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 Å−1 employing the vV algorithm reported 2.78 × 10−9 m2/s obtained from MSD (mean square displacement, equation (9)) and 2.75 × 10−9 m2/s from VAC (velocity autocorrelation, equation (10)) function at 300 K, which overestimate the experimental data at the same temperature (2.39 × 10−9 m2/s [26, 27] and 2.49 × 10−9 m2/s [26–28]). However, TIP4P/2005 water model [29] gives an excellent agreement, 2.39 × 10−9 m2/s, obtained from both the MSD and VAC in our previous study [25] with the experiment one.Table 6
Water diffusion coefficients (10−9 m2/s) at 25°C employing a velocity Verlet algorithm, where n is the number of (Na+ + Cl−) ions.
System (n)κ = 0.05 Å−1κ = 0.10 Å−1κ = 0.25 Å−1A1 (2)2.742.812.68A2 (4)2.762.792.66A3 (8)2.802.822.64A4 (16)2.872.822.62A5 (36)3.012.842.52A6 (72)3.152.802.38A7 (144)3.252.582.05A8 (216)3.142.321.76A9 (288)2.942.021.51Dw 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 Dw employing the vV algorithm underestimates except for ST and Ewald with κ > 0.45 Å−1. These MD simulations [8] reported that Dw is 3.03 × 10−9 m2/s in the Na+-water system and 3.09 × 10−9 m2/s in the Cl−-water system using ST employing the Gear algorithm with Nw = 215 of SPC/E water. An MD simulation for alkali Earth metal cations (Mg2+, Ca2+, Sr2+, and Ba2+) in an aqueous solution at 25°C of SPC/E water potential (Nw = 215) using Ewald employing the Gear algorithm reported that Dw is in the range of 2.46–2.62 × 10−9 m2/s [17].In Table6, roughly speaking, Dw with κ = 0.05 Å−1 increases up to the A7 MD system with slow decreases as c increases, while Dw with κ = 0.10 Å−1 has almost the same values up to the A6 MD system with sudden decreases, and Dw with κ = 0.25 Å−1 decreases monotonically up to the A6 MD system with sudden decreases. The behavior of Dw with κ = 0.10 Å−1 seems to be the most reasonable one: Dw 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 D0 = DPBC + 2.837297 kBT/(6πηL), where DPBC is the diffusion coefficient calculated in the simulation, kB 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, D0-DPBC = 0.1856 × 10−9 m2/s which is quite small, compared to our calculated Dw = 2.02–2.84 × 10−9 m2/s with κ = 0.10 Å−1 (Table 6). This correction for Dw 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, Dw, decreases with increasing salt concentration and for other structure-breaking salts, Dw increases with increasing salt concentration c. Both ratios of experimental and simulation Dw (NaCl)/D0 (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 Å−1 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 Å−1 employing a velocity Verlet algorithm, where n is the number of (Na+ or Cl−) ion(s).
System (n)Na+-Na+Cl−-Cl-Na+-Cl-Real reciprocalReal reciprocalReal reciprocalA1 (1)————−0.293.83A2 (2)0.0013.40.9315.7−11.0−2.16A3 (4)1.7030.92.8926.5−23.3−7.59A4 (8)6.8633.26.6734.9−48.2−18.0A5 (18)34.138.029.939.3−100−42.8A6 (36)52.834.542.731.7−131−32.9A7 (72)10715.990.214.1−193−13.5A8 (108)18214.416212.8−272−11.9A9 (144)2574.912324.61−347−4.05The 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. Figure8 shows the sum of diffusion coefficients (○) of all ions multiplied by n in DA1 to DA9 in Table 3 for systems A1 to A9 for κ = 0.10 Å−1. The sum of the “Ideal” diffusion coefficients (---) of all ions obtained from DA1 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 Table3 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 Λm 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 gii (r) with κ = 0.25 Å−1 show the ordinary behavior similar to the ion-water radial distribution functions giw (r). However, those functions gii(r) with κ = 0.10 Å−1 give a dramatic change with a very deep minimum of gNaCl(r) and, as a result, sharp maxima of gNaNa (r) and gClCl (r) at the distance R1 = 9.95 Å. The ratios of nNaCl (R1)/(nNaNa (R1) + 1) and nClNa (R1)/(nClCl (R1)+1) for the A5 ∼ A9 MD systems with κ = 0.10 Å−1 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 Å−1 gives the closest one to the experimental data. The behavior of Λm 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 ofDNa+ and DCl− for NNaCl = 1, 2, 4, 8, and 18 are obtained with κ = 0.25 Å−1. 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 Å−1 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 DNa+ and DCl−. We conclude that the appropriate value of κ is determined by the number of ions rather than the water molecule.



---

*Source: 6661196-2020-12-09.xml*