Effects of ion concentration and solvent composition on the properties of water-methanol solutions of NaCl. NPT molecular dynamics computer simulation results

Isothermal-isobaric molecular dynamics simulations are used to examine the microscopic structure and other properties of a model solution consisting of NaCl salt dissolved in water-methanol mixture. The SPC/E water model and the united atom model for methanol are combined with the force field for ions by Dang [J. Amer. Chem. Soc., 1995, 117, 6954] to describe the entire system. Our principal focus is to study the effects of two variables, namely, the solvent composition and ion concentrations on the solution's density, on the structural properties, self-diffusion coefficients of the species and the dielectric constant. Moreover, we performed a detailed analysis of the first coordination numbers of the species. Trends of the behaviour of the average number of hydrogen bonds between solvent molecules are evaluated.


Introduction
Theoretical studies of ionic solutes in water have generated an enormous amount of literature during the recent decades. However, several problems still require a comprehensive solution, particularly due to the lack of a complete understanding of water as a solvent. The recent progress in the area has been related to the application of computer simulation methods to the models of various charged species in aqueous media with different degrees of sophistication, see, e.g. [1,2].
Less studies have been devoted to the properties of ionic solutes, even simple salts, in combined solvents that include water and an organic component. These systems exhibit a wide range of physicochemical phenomena. They are of great importance for basic and applied research, but not easy to deal with due to their complexity. The amount of experimental data in this area is much less comprehensive compared to ionic aqueous solutions.
For specific purposes of our project, it is worth mentioning that in a large number of works the molecular dynamics computer simulations have been applied to mixtures of water with various solvents.
The mixtures of water with alcohols, with DMSO, DMF, acetonitrile, acetone and glymes seem to be most frequently studied, see, e.g. [3][4][5][6][7][8] and references therein. Actually, practical needs, besides academic interest, determine the choice of systems for experimental research and for computer simulations.
The principal objective of the present work, as a part of an ampler project, is to investigate the microscopic structure and other properties of sodium chloride ionic solutions at different ion concentrations with water-methanol mixed solvent in the entire range of compositions. In spite of impressive efforts and success of long-term studies coming from the laboratory of Hawlicka [11][12][13][14][15][16][17][18][19], there are several issues that require further exploration of this type of solutions. Some of these issues will be discussed below.
Our study is performed in the framework of molecular dynamics computer simulations. The structure of solutions is given in terms of different descriptors, namely, a brief account of the radial distribution functions is provided, and the first coordination numbers are discussed. In addition, we would like to analyze the conditions under which the cluster formation of ions occurs. With this aim, we analyze the time evolution of the number of clusters and other related distributions. The structure of the maximal cluster formed at a certain ion concentration and solvent chemical composition is visualized. Insights into the dynamic properties are obtained by exploring the mean square displacements as a function of time and the resulting self-diffusion coefficients of the species. Changes of hydrogen bonds statistics are analyzed. Finally, we discuss the behaviour of the dielectric constant of the solutions upon the changes in the solvent composition and ion concentration.
At this initial stage of the project, we performed all the analyses using a single combination of the force field models for water, methanol and ions as well as a single combination rule for the cross interactions. Sodium chloride, well soluble in water, has been chosen here as a test case to investigate other salts and co-solvents in our future work. Moreover, our study is restricted to room temperature and ambient pressure (1 bar). Only the changes of salt concentration and the chemical composition of the solvent are investigated. Wider, more precise and profound insights into the thermodynamics of this type of solutions, the analyses of an ampler set of dynamic and dielectric properties, require additional explorations.

Model and simulation details
Prior to discussing the simulation techniques, we would like to comment on the modelling used in the previous studies of alkali halides solutions with water-methanol solvent. The system in question, i.e., NaCl-water-methanol, is complicated for simulations as it involves a model for each compound. Moreover, certain combination rules should be used for cross interactions.
In almost all cases, Hawlicka et al. [11][12][13][14][15][16][17][18][19] used the three site Bopp-Jancso-Heinzinger (BJH) model for water [28] and a three-site Palinkas-Hawlicka-Heinzinger (PHH) model for methanol [29]. The ionsolvent and ion-ion potentials were derived from ab initio calculations by fitting the energies to analytical expressions, see, e.g. [30]. In some studies, the Lorentz-Berthelot (LB) combination rule was used. All the simulation results published by Hawlicka group have been obtained using either microcanonical (NVE) or canonical (NVT) ensemble, the simulation box sizes in each procedure have been adjusted to the experimental density values, that unfortunately have not been specified but in a few cases turned out to be in terms of simulation box length.
Recent report by Bouazizi and Nasr [26] combined the flexible SPC model for water and the flexible six-site model for methanol (apparently, the OPLS all-atom model) with Smith and Dang model for ions [31]. The LB combination rules were used. The results of the procedure follow from the NVT simulations using Berendsen control for temperature, where a rather small number of particles (three ion pairs and 250 solvent molecules) was simulated. A few runs were performed with 1024 solvent molecules.
The most recent study of the system in question performed in [32] involve the SPC/E model for water [33] as well as the models for methanol and for ions developed in the laboratory of Vrabec [34,35]. Again, the LB combination rules were employed. However, a correction into the cross-diameter rule was adopted for water-methanol. The simulations were performed in the isothermal-isobaric (NPT) ensemble.
Our calculations have been performed in the isothermal-isobaric (NPT) ensemble at 1 bar, and at a temperature of 298.15 K. We used the GROMACS software package [36], version 5.1.2. For methanol, we used the united atom model [37]. The SPC/E model was taken for water [33]. Lennard-Jones parameters for sodium and chloride ions were chosen according to the force field [38]. The Lorentz-Berthelot combination rules were used to determine cross parameters. The nonbonded interactions were cut off at 1.4 nm. The long-range electrostatic interactions were handled using the particle mesh Ewald method implemented in the GROMACS software package (fourth order, Fourier spacing equal to 0.12) with a precision of 10 −5 . The van der Waals tail correction terms to the energy and pressure were taken into account. In order to maintain the geometry of the water and methanol molecules, the LINCS algorithm was used.
As concerns the procedure, a periodic cubic simulation box was set up for each system. The GRO-MACS genbox tool was employed to randomly place all particles in the simulation box. To remove possible overlaps of particles introduced by the procedure of preparation of the initial configuration, each system underwent an energy minimization using the steepest descent algorithm implemented in the GROMACS package. Minimization was followed by a 50 ps NPT equilibration run at 298.15 K and 1 bar using a timestep of 0.25 fs. We used the Berendsen thermostat and barostat with τ T = 1 ps and τ P = 1 ps during equilibration. Constant value of 4.5 · 10 −5 bar −1 for the compressibility of the mixtures was employed. In the case of pure methanol solvent, the compressibility was taken to be 1.2 · 10 −4 bar −1 . The V-rescale thermostat and Parrinello-Rahman barostat with τ T = 0.5 ps and τ P = 2.0 ps and the time of step 2 fs were used during the production runs. Statistics for each mole solvent composition and various ions concentration for any of the properties were collected over several 10 ns NPT runs, each starting from the last configuration of the preceding run. The time extension for each series of calculations will be mentioned below in the appropriate place but not less than 70 ns.
In order to compare the density of the solutions with experimental results given on the molality scale, see the subsection 3.1 below, the total number of particles (ions and molecules) in the box ranged from ≈ 2500 to ≈ 3700, the number of NaCl ion pairs was in the interval from 0 to 40.
On the other hand, while exploring the composition changes in the solvent, the total number of solvent molecules was kept fixed at 3000. Three series of calculations correspond to "inserting" the N α (α = Na, Cl) of positive and negative ions (N α = 20, 40 and 60) into the box containing 3000 solvent molecules. Hence, the ion fraction, X ion , [X ion = 2N α /(N m + N w + 2N α )] is equal to X ion ≈ 0.013, X ion ≈ 0.038 and X ion ≈ 0.0625, respectively. The composition of the solvent was being changed in each series of calculations; it is described by the mole fraction of methanol molecules X m [X m = N m /(N m + N w )]. In the case X m = 0, i.e., in a pure water solvent, the prepared ionic solutions correspond to the concentration 0.366 mol/dm 3 , 1.081 mol/dm 3 and 1.77 mol/dm 3 , respectively. Additional series of calculations have been performed for the salt-free case with 3000 molecules at different X m . In this case, in the absence of ions, it was sufficient to terminate the runs at the time of 60 ns.

Density of solutions
We begin with the description of the effect of the density of solutions on the solvent composition and on the ion concentration. The simulation results are compared with the experimental data for sodium chloride aqueous solutions and with the data reported by Takenaka and Takemura [39] for solutions involving a solvent of two species, water and methanol. All the data are collected in table 1 below.
The principal conclusion made from the data shown in the table is that the model used in the present study leads to the results that favourably agree with experimental data of [39] in a wide interval of molality at various compositions, namely, at X m = 0.129; 0.307; 0.571 and 0.917. At a high fraction of methanol as a co-solvent, the discrepancy between the simulation results and experimental data is more pronounced. However, it is of interest to perform similar, but more extensive, comparisons for two limiting cases of a solvent composition, namely, for pure water as a solvent (X m = 0) and for pure methanol as a solvent (X m = 1), since the experimental data are available. In the case of pure water as a solvent, we performed simulations for NaCl solutions at a molar salt concentration reported in table 2 of [41] using a much larger number of molecules in the box, as mentioned in the previous section. Our simulation data are given in the first column of the above table. The experimental values for these conditions (converted to molality scale) were derived by interpolation from the data presented by Romankiw et al. [40]. The reason for doing that is due to the fourth line in A further confirmation of this behaviour is illustrated in the left-hand panel in figure 1. Both the model of the present study and the model used by Kohns et al. [32] yield a very good agreement with the experimental results of [43,44] at low and intermediate molalities. At high values of molality, our model underestimates the density of the solution whereas the model by Kohns et al. [32] overestimates the density. Finally, the model employed in [41] is in agreement with the experimental data of [39] at low molality and works reasonably well at intermediate values of molality. In the right-hand panel in figure 1, we use the experimental data from [43,45] for the methanolic solutions of NaCl. Our simulations correctly reproduce the changes of the density of these solutions depending on molality. The model in question slightly overestimates the density but the discrepancy between simulations and experimental results is of the order of a fraction of one percent of the value for the density. Having this satisfactory description of the solutions in two limiting cases, we have undertaken simulations for NaCl in mixed solvents, figure 2. Here, we show the results for the salt-free case and for three sets of systems with a varying solvent chemical composition. Each set corresponds to a different fixed number of ion pairs, N ip = N Na = N Cl , added to the solvent, namely, with N ip = 20, N ip = 40 and N ip = 60. In addition, we reproduced a piece of the curve for a set of systems with N ip = 50. In all cases, the number of solvent molecules is fixed at 3000, i.e., N w + N m = 3000. The composition dependence of the density of the solution for systems with N ip = 20 is similar to the set of salt-free systems for all X m . Interestingly, the 23601-4 exp. data Ref. [39] simulation Ref. [45] simulation Ref. [40] present work NaCl in water solvent exp. data Ref. [43] [32] and from [41], respectively. The experimental results marked by circles are from [39]. The experimental data marked by diamonds are from [44]. Panel (b): dependence of the density of NaCl solutions in methanol solvent (united atom model) on ion concentration from simulations of the present work. Experimental data are from [45] and [43], circles and squares, respectively. "insertion" of this amount of ion pairs into a mixed solvent results in a uniformly augmenting density for all X m , including two limiting cases, X m = 0 and X m = 1. A further increase of the amount of cations and anions in the system (N ip = 40, N ip = 50 and N ip = 60) leads to higher values of the density of the solution for various compositions in terms of X m . However, each of these three lines has a terminating point at which the systems reach saturation. Three stars in figure 2, if joined together, describe a part of the saturation curve. The star upon each line describes the "last" homogeneous system prior to the formation of a single rather big cluster (involving the majority of ions, which permits to consider it as a nucleus or a grain of a solid crystal phase) at a higher value of X m . This cluster coexists with a solution, i.e., with the entire solvent subsystem containing much less ions that are either solvated cations and anions or ion pairs. We will return to this issue in the body of the manuscript below.

Pair distribution functions, coordination numbers and hydrogen bonding analysis
The microscopic structure of the solutions under study is described in terms of various pair distribution functions (pdfs). Many of them were described in every detail in several publications for the aqueous solutions of NaCl (X m = 0), see, e.g. [41], as well as for NaCl solutions with a composite water-methanol solvent, see, e.g. [12][13][14][15]26]. To avoid unnecessary repetitions, we refer to the most comprehensive and to the most recent description of the pdfs given in [26]. Moreover, the discussion by these authors involves comparisons with previous developments and with a few experimental findings concerning a microscopic structure of salt-free mixtures and some solutions.
Our few examples, shown in figures 3 and 4, just serve to confirm that the model of the present study leads to correct and physically well grounded conclusions in accordance with previous developments. The results are given solely for a single set of systems, namely the number of solutes is fixed at N ip = 40 while the solvent composition changes with X m , the total number of solvent molecules is kept fixed at N w + N m = 3000. As concerns the behaviour of the cross pdfs between solvent species, g OW-OM (r) in figure 3, we observe that the first maximum gradually grows whereas the second maximum monotonously decreases with an increasing amount of methanol in the solvent. Thus, the water-methanol contacts become more probable with an increasing X m . However, the correlations between two species at a larger scale become weaker. This behaviour is in accordance with the trends observed for the salt-free systems. At the same time, the short-range correlations within a water subsystem become stronger with an increasing X m , as evidenced by the growth of the first maximum of g OW-OW (r) shown in the right-hand panel in figure 3. The positions of the first maxima describing most probable configurations of solvent species in terms of g OW-OM (r) and g OW-OW (r) are practically unchanged with the changes of X m . Most probably, the correlations between more distant molecules are affected by changes of the solvent chemical composition, see [26] for an ampler discussion.
As concerns the behaviour of ion-water pdfs, we would like to mention the following. In the water-rich solvents, see, e.g. the curves corresponding to g ion-OW (r) at X m = 0.1, ions are very well distributed in the solvent, figure 4 (left-hand panel). This is evidenced by a high first maximum of g Na-OW (r) and a well pronounced second maximum of this function. On the other hand, the height of the first maximum of g Cl-OW (r) and of g OW-OW (r) are nearly equal. Inspection of the curves describing g Na-OM (r) and g Cl-OM (r) (right-hand panel in figure 4) and their comparison with the g Na-OW (r) and g Cl-OW (r) functions leads to the following conclusions. The methanol molecules are capable of approaching the sodium cation as satisfactorily as the water molecules do, in spite of the number of methanol particles being quite low at X m = 0.1. On the other hand, the chloride anion is preferentially surrounded by waters, the methanol species (OM) possibly enter the solvation shell of Cl − , though with much lower probability compared to OW. In solvents that are rich in methanol (X m = 0.9), we observe that sodium cation attracts as much waters as available, but the solvation shell contains some methanol molecules since the first maximum of g Na-OM (r) is quite high. The chloride anions also attract as many water molecules as available at this solvent composition, but the magnitude of the first maximum of the function g Cl-OM (r) substantially grows upon the changes of X m from 0.1 to 0.9. The ion-ion pair distribution function, g Na-Cl (r) (in figure 5) exhibits a maximum at r ≈ 2.8 Å describing the probability of finding contact ion pairs, the second maximum at r ≈ 5 Å showing the presence of solvent-shared ion pairs, while the following oscillations are not visible on the scale of the figure. The first and the second maxima essentially increase in magnitude with the change of the solvent composition. At the same time, the pdfs describing the configurations of similarly charged ions are characterized by a more pronounced structure [we show the function g Na-Na (r) at X m = 0.9]. It is worth mentioning that in the case of dilute solutions, the trends of the behaviour of the ion-ion pdfs with ion concentration are in accordance with the calculations of the potential of the mean force, see, e.g. [46]. In the present case, at X m = 0.9, we observe the development of anti-phase oscillations that witness the formation of short-range ordering of ions, which results in cluster formation upon a further increase of the methanol content.
We should like to finish the discussion of the pair distribution functions with the following comment. We believe that the model yields a qualitatively correct and a physically sound picture, while other combinations of the force fields yield formally similar results. However, any kind of improvement could be attempted only if the experimental structure factors of such complex systems were available. Then, the fitting of computer simulation results to, e.g., diffraction experiment data, in spite of the expected technical difficulties, could be performed along the lines proposed in [47,48] for aqueous electrolyte solutions.
The pair distribution functions yield the running coordination numbers through the equation, where ρ j is the number density of species j. The first coordination number is obtained by putting R = r min , i.e., at the first minimum of the corresponding pair distribution function. As concerns the coordination numbers describing the solvent species, we can observe that n OW-OW coordination number almost monotonously decreases if X m grows from zero to unity, figure 6. It starts at a value ≈ 4.5 close to what is expected for water. At the same time, n OW-OM increases. The corresponding line ends up at ≈ 3.2, indicating a less coordinated water molecule by methanol particles in comparison with the aqueous medium. We were unable to detect pronounced overall differences between the coordination numbers of solvent species upon changes of the number of solute ions from N ip = 20 to N ip = 60. Just in the latter case the lines deviate from the results corresponding to N ip = 20 and approach the saturation conditions at X m ≈ 0.6. The shape of the lines joining the calculated points is approximately linear.  The behaviour of the coordination numbers of ions and solvent species exhibits certain interesting features and peculiarities. Namely, from the left-hand panel in figure 7 (N ip = 20), we learn that n Na-OW linearly decreases with increasing X m whereas the coordination number n Na-OM almost linearly increases. The coordination number n Cl-OW decreases as well, upon increasing X m . However, it follows the trends of n Na-OW until X m ≈ 0.4. In the interval X m > 0.6, the decay of n Cl-OW is much faster. A qualitatively similar shape with two different rates of changes is seen in n Cl-OM . The difference in coordination of Na and Cl at X m = 0 is in accordance with the previous studies of aqueous solutions. However, we can observe that the difference of coordination of sodium cation in pure water and in pure methanol is not big. On the other hand, the chloride anion is much less coordinated in methanol compared to water.
The values for the coordination numbers at X m = 1 in the left-hand panel in figure 7 agree with the results obtained recently in [49], see table 9 of that reference. Similar trends of behaviour of ionsolvent species coordination numbers are observed for a higher number of solute ions, N ip = 40. While approaching saturation conditions, the coordination numbers of cations and anions in the case N ip = 40 become close to each other, resembling what happens in the systems with a lower number of solutes in the solvent with a predominating methanol amount. However, the system with N ip = 60 can be very close to saturation, but the coordination numbers of sodium cations and of chloride anion differ much (see dotted curves with squares in the right-hand panel in figure 7). Thus, the ion-solvent coordination numbers do not indicate on their own that the system approaches saturation.
Finally, the first coordination number, n Na-Cl , as a function of the solvent composition is given in figure 8 for three sets of systems differing by the ion fraction. All the curves behave qualitatively similarly: the coordination number grows with increasing X m and exhibits a substantial jump at conditions when a rather big ionic cluster is formed. The magnitude of the jump depends on the number of solute molecules as indicated in figure 8. Thus, the n Na-Cl exhibits a single and well pronounced peculiarity in a water-methanol solvent.
The threshold value for the observed jump depends on the solvent composition (one may express it as X m or in terms of weight percentage of the co-solvent) and on the ion concentration (usually molarity or molality scale is used in experimental studies). Certainly, the precise values for threshold conditions depend on the details of the force field of ions as well as on the type of the models for water and methanol, e.g., united atom versus all-atom modelling. Therefore, we postpone a more detailed exploration of this issue in terms of the coordination numbers to a future work. However, it is worth mentioning that the threshold situation can be monitored by the behaviour of the self-diffusion coefficients of ions and by the conductivity of the solutions, besides the coordination numbers.
Our final concern in this section is in describing the changes of hydrogen bonds in a solvent subsystem  of the solutions under study upon composition changes. The average number of hydrogen bonds per water molecule, as a function of X m , is shown in figure 9. Technically this was obtained using the gmx hbond utility with default options of the GROMACS software.
The curves are rather smooth and exhibit certain similarities to the behaviour of coordination numbers discussed above. Furthermore, for the sake of comparison, we performed necessary calculations for the salt-free systems of variable composition. The limiting values for water (X m = 0) and for methanol (X m = 1) are in accord with the previous study from this laboratory [4].
In the entire interval of X m , the average number of H-bonds per water molecule n HB is lower in the system with ions in comparison with a salt-free system at the same X m . The cumulative effect of solutes (cations and anions) is that they actually make the formation of a certain amount of water-water H-bonds impossible due to the existence of the solvation shells. By considering the magnitude of changes, this cumulative effect can be referred to as primary effect. At a fixed X m , the average numbers of the cross water-methanol H-bonds (per the number of water molecules) in the combined solvent, i.e., in the absence of ions, and in the solution are close to each other. Finally, even a weaker effect is observed on the average number of bonds formed between methanol molecules. From the analyses above, it follows that ions affect the properties of the system mainly via the modification of the structure of a water subsystem of the solvent. As a result, the trends of mixing of water with methanol may display certain changes but not big ones.
Final part of this subsection is devoted to the description of conditions under which the ion clusters are formed. This problem has already been considered for NaCl aqueous solutions in several publications, see, e.g., [50][51][52]. The cluster (or clusters) formed during the final part of a long trajectory obtained in simulations is the precursor or a small grain of the salt precipitate. The study of the solubility of a given salt in a solvent of different composition is beyond the scope of the present work at this stage of our project. Nevertheless, we have picked up the final frame of the MD trajectory for the system with N ip = 60, X m = 0.6333 (N w + N m = 3000) and made visualization of the coordinates of particles. It is given in figure 10 (water molecules in the box are shown besides the ions). Apparently, a nice picture with a pronounced symmetry of ions in a rather big cluster does not provide a description of the structure formed in quantitative terms [53]. Still, one can appreciate that not all the ions belong to the cluster. On the other hand, in spite of macroscopic miscibility, heterogeneity of the distribution of solvent species on the local scale (say, the scale of the simulation cell) is seen, methanol molecules fill the "empty" space where water particles and ions are not present.
In order to get a better and more quantitative insight into the formation of a big cluster of ions, we plot the number of clusters, the maximal cluster size and the average cluster size as functions of time emerging from the trajectory of the system. Technically, these properties result from the gmx cluster utility. The parameter for distance to count the ions belonging to a cluster is taken a bit larger than the location of the first maximum of the g Na-Cl (r). A bit different choice can yield slightly different numbers, but qualitatively the trends in the behaviour of the properties describing clusters remain the same.
Inspecting figure 11, we can observe that prior to saturation conditions, i.e., at X m = 0.6, the number of clusters fluctuate in the interval approximately between 100 and 110 (recall that we have 120 ions). In other words, the majority of ions are "free" and only some of them presumably form ion pairs along the entire trajectory. On the other hand, crossing the saturation leads to drastic changes of this function. After approximately 20 ns, the number of clusters drops to reach a stable value around 30 clusters. The maximal cluster size, in terms of the number of ions participating in it, becomes around 90, the average cluster size attains a much bigger value for X m = 0.633 in comparison with X m = 0.6. Thus, the majority of ions can really be considered as those belonging to a rather big cluster, the rest of the ions can be free or be involved in small entities such as ion pairs, for example.

Self-diffusion coefficients and the dielectric constant
The self-diffusion coefficients of ions, i.e., water and methanol, in our work were calculated from the mean-square displacement (MSD) of a particle via Einstein relation, where i refers to water or methanol, or ions Na, Cl, and τ is the time origin. Default settings of GROMACS were used for the separation of the time origins. Moreover, the fitting interval (from 10% to 50% of the analyzed trajectory) has been used to calculate D m and D w . Also, a special care has been taken regarding the fitting for the cases with a small number of particles on the extremes as well as of cations and anions along the X m axis. A set of our results is given in figure 12. First, it is worth mentioning that we performed additional calculations for the salt-free system in order to establish the magnitude of the effects due to the presence of ionic solutes. The D w is lower for all systems with ions (X ion 0) compared to the salt-free system (X ion = 0) in the entire solvent composition interval, upper left-hand panel in figure 12. The result for D w , in a salt-free system at X m = 0, is undoubtedly correct, which is around 2.6 (for SPC/E model). A bigger amount of ion solutes leads to a reduction of the value for D w . However, a weaker effect of this sort is observed in water-rich solvents, apparently because there is enough water molecules to provide solvation of all the ions in the system and still there is enough "free" water. A more pronounced reduction of D w is observed in the systems at intermediate values of X m , because practically all water molecules belong to the solvation shell of ions.
As concerns the self-diffusion coefficient of methanol species (right-hand upper panel in figure 12) we would like to mention the following trends. For the salt-free system, the value for D m is correct, which is around 2.1 at X m = 1. The ions are solvated by methanol molecules in the methanol abundant region of the solvent compositions. Consequently, D m is lower for all such solutions in comparison with the salt-free system. If the number of solutes is low (N ip = 20), the self-diffusion coefficient D m in the solution is practically equal to the values obtained for the salt-free system in case X m is low (in the water-rich solutions). In other words methanol molecules do not permeate the solvation shells of ions at these conditions. Nevertheless, a reduction of the values for D m is pronounced if the amount of solute ions increases from N ip = 20 to N ip = 40 and 60. In particular, if the number of methanol molecules is low, i.e., at low values of X m , the reduction of D m is due to hydrogen bonds between methanols with water molecules in the solvation shells as well as due to methanols "directly" participating in the solvation of ions. It is impossible to precisely evaluate the contributions from each of the two factors separately.
An overall picture of the dependence of self-diffusion coefficients of all species in the solution on the composition is shown in the lower left-hand panel in figure 12. The behaviour of D w and D m on composition resembles the trends already observed for the salt-free mixtures with a minimum of each function in the region of intermediate X m , see, e.g. [4,13]. On the other hand, we are able to definitely state that D Cl and D Na are both much lower than D w and D m . Moreover, we are convinced that D Cl is higher than D Na in the entire interval of X m , and, in particular, in the limit of pure water as a solvent. The results obtained after the runs of 90 ns still do not provide quite accurate values. In spite of the error margins, we note that both D Cl and D Na have a minimum on a methanol-rich side of the solvent composition. To our best knowledge, an alternative method for calculation of D i via the velocity autocorrelation functions provides the same accuracy. It is surprising that the curves obtained by this method with only three ion pairs (and with 250 solvent molecules) for the system in question, but with the use of a different model, are nearly perfect (see figure 8 of [26]). Still another modelling employed by Hawlicka et al. [13] with solely four data points for solutions with a mixed solvent indicates that both D Cl and D Na have a minimum at X m = 0.5. Therefore, it is clear that further studies of the systems with a larger number of particles are indispensable.
As concerns the effect of ion concentration, we attempted to evaluate the self-diffusion coefficients with a larger number of ions, N ip = 40, and to make comparisons with the previous case, N ip = 20. These results are shown in the right-hand lower panel in figure 12. Both, the D Cl and D Na , decrease in magnitude at a higher ion fraction. The reason of such a behaviour is in a higher density of well-solvated ions in the system. Consequently, their movement is more hindered.
The last part of this subsection is concerned with a few peculiar issues of the dielectric response of NaCl solutions with water-methanol mixed solvents. Since an electrolyte solution is a conductor, the mobile ions move in response to an imposed static electric field. Thus, while dealing with this type of systems, it is common to define a static equilibrium property that plays a role similar to the dielectric constant of an insulating fluid, see, e.g., a comprehensive discussion of the problem in [54]. Actually, the apparent dielectric constant reported from the experimental measurements is an extrapolation from the frequency dependent dielectric constant, as documented, for example, in [55] while discussing the dielectric properties of NaCl aqueous solutions of different salinity. However, one should have in mind that there exist additional dynamic contributions hidden in the extrapolated values. In several publications that describe the application of molecular dynamics computer simulations to electrolyte solutions, the static dielectric constant is defined or estimated in terms of the time-average of the fluctuations of the total dipole moment of the system [56][57][58][59], obtained by the summation of the dipole moment vectors of the entire set of water and methanol molecules in the simulation box (see, e.g., the GROMACS manual available online and [36]), where k B is Boltzmann's constant and V is the simulation cell volume. The average value of all Cartesian projections of M should be small to provide the evidence that particles on the scale of the box are distributed homogeneously, see, for example, figure 5 (c) of [60] for illustration of similar calculations for pure water. We use the procedure based on equation (3) to get an insight into the predictions of the model that involves a specific set of force fields (for ions, water and methanol).
As discussed in detail in [61] for aqueous electrolyte solutions with completely dissociated solute molecules (NaCl in water is an example, see table 7.6 of [61]), the static permittivity of the sample, ε, decreases with increasing the electrolyte concentration. This permittivity simultaneously represents the static permittivity of the solvent in the solution. Moreover, the concentration dependent decrease of the static solvent permittivity is defined as the dielectric depression and is empirically well fitted by the linear dependence. In the present case, the situation is more complex since the dielectric depression is determined by the solvent composition, apart from the ion concentration. Furthermore, at the moment we are unaware how sensitive the dielectric depression is to the details of the force fields for water, methanol and ions. The curves from simulations for a combined water-methanol solvent and for the three solutions at different ion fractions are shown in figure 13. As it follows from the comparisons of the dielectric constant of a salt-free system and the experimental data performed in the recent work from this laboratory [4], the solvent model underestimates the values for ε in the entire interval of the composition. Nevertheless, the curve behaves qualitatively correctly, and the discrepancies between the simulation results and experimental data are not big, see figure 4 of [4]. Principal effects of ions can be seen on the waterrich side. Namely, the dielectric constant decreases in magnitude with an increasing ion fraction. This behaviour is qualitatively correct. On the other hand, on the methanol-rich side, all three curves describing the systems with different ion fractions flow closer, indicating a weaker response of the systems with a predominant number of methanol molecules. Undoubtedly, it is of interest to extend the study of dielectric properties along various lines briefly discussed above in a specifically focused future research.

Summary and conclusions
We have performed an extensive set of molecular dynamic simulations in the NPT ensemble to study properties of NaCl solutions with water-methanol solvent in the entire range of a solvent composition. We explored systems at a fixed number of solvent molecules and three values of the number of ion solutes. The salt-free system has been discussed as well. Out of many possible versions of the model for a solution, only one was studied. We consider this as a first step of more systematic studies of solutions involving ionic solutes of different complexity and a solvent mixture of water and methanol at different thermodynamic states using nonpolarizable models.
We explored the evolution of a microscopic structure of all the systems in terms of the pair distribution functions together with the first coordination numbers of the species. We approached saturation at certain thermodynamic conditions. The number of clusters and maximal cluster size as a function of time have been explored. Statistical features of the hydrogen bonds in a water subsystem and of water-organic co-solvent bonds have been described. The self-diffusion coefficients of all the species and the dielectric constant were also calculated. All the simulations were performed at room temperature and at ambient pressure, 1 bar.
From a comparison with a limited set of experimental data for the systems in question and with the results of other authors on the related systems, we can conclude that the predictions obtained are qualitatively correct and give a physically sound picture of the properties explored. We observed that practically all the properties investigated are sensitive to the composition of a water-organic liquid solvent.
At the present stage of the development in the field, the missing and very important elements are really numerous and need much more detailed investigations. Namely, it is necessary to investigate other force fields for methanol species at the united atom and/or all atom level of modelling for the type of solutions of this work. A comprehensive discussion of the modelling strategies that focused onto a better understanding of aqueous solutions has been recently presented by Nezbeda et al. [1]. The role of combination rules should be explored as well, see, e.g., [51] for the case of aqueous solutions of NaCl. A possibility for fitting parameters of different models to avoid a spontaneous cluster formation in accordance with experimental findings for stable solutions with composite a water-methanol solvent has not been explored so far, in contrast to simpler aqueous solutions [62]. It is extremely desirable to perform comparisons with structure factors coming from experiments (we are not aware of the experimental data for NaCL solutions with water-methanol solvents) in close similarity to the developments available for salt-free mixtures, e.g., [3], and, e.g., for aqueous solutions of salts with a different complexity [63][64][65][66]. On the other hand, trends of the behaviour of thermodynamic, and, in particular, mixing properties require further investigations. We hope to address some of these issues in our future studies. The present work has partly paved the way toward the understanding of what should be done.