On the properties of methanolic NaCl solution by molecular dynamics simulations

Isothermal-isobaric molecular dynamics simulations are used to examine the microscopic structure and principal thermodynamic properties of a model solution consisting of NaCl salt dissolved in methanol solvent. Four united atom force fields for methanol are involved. Concerning ion solutes we used the Joung-Cheatham, Smith-Dang models as well as the model from the laboratory of Vrabec. Our principal focus is to evaluate the quality of predictions of different combinations of models for basic properties of these solutions. Specifically, we explored the change of density on molality, the structural properties in terms of various pair distribution functions, the coordination numbers, the number of ion pairs and the average number of hydrogen bonds. In addition, changes of the self-diffusion coefficients of species, the solvent dielectric constant and the evolution of the surface tension with ion concentration are described.


Introduction
From a general perspective, this work is focused on the theoretical description of one simplest example of a class of electrolyte systems that comprise ionic solutes dissolved in water with organic co-solvent. Molecular dynamics computer simulations are used as tools to explore the microscopic structure, ions solvation and thermodynamic properties of interest.
The present contribution is a part of our ongoing project focused on the theoretical exploration of properties of ionic solutes in water-alcohol solvents of variable composition. At the first stage, watermethanol model mixtures have been explored in this laboratory [21][22][23][24], continuing several previous studies of this kind of systems, see e.g., [25][26][27][28][29][30][31][32] and respective references in the above cited works. Still, there remains room for additional insights, having in mind the recent novel modelling of liquid methanol [33,34].
One limiting case of ionic solutes in water without organic co-solvent generated huge amount of computer simulation results during past decades. Several remaining problems within this methodology are due to a limited understanding of water as solvent on its own, and of the influence of ions on water behaviour. For example, modelling of water in the framework of a particular model makes necessary to use adequate, specific set of parameters describing force field for ionic solutes. This issue is well documented for example in the development of frequently used Joung-Cheatham (JC) [35] model for alkali halides aqueous solutions with TIP3P [36], TIP4P-Ew [37] and SPC/E [38] water models. Methodological aspects of the development and application of ionic force fields to NaCl aqueous electrolytes have been critically reviewed in [39][40][41][42]. The state of art for this particular aqueous solution is well described in [43][44][45].
On the other hand, the most comprehensive investigation of ionic solutions in water-methanol solvents was performed in the laboratory of Hawlicka [4][5][6][7][8][9][10][11][12]. These studies were restricted to NaCl and NaI univalent salts, as well as MgCl 2 and CaCl 2 salts with divalent cations as solutes. The solvent species was described by using the BJH water [46] and PHH [47] methanol models. In addition, the ab initio calculations for ion-solvent interaction energies were involved. Some recent contributions contributed to a better knowledge of the properties of NaCl solutions in water-methanol solvent as well [48,49].
In contrast to aqueous NaCl solution, the non-aqueous NaCl solution with pure methanol as solvent is much less studied. Even the amount of experimental data concerning the properties of this system is much less comprehensive. A set of models for liquid methanol, designed for application within the computer simulation methods, is quite ample [33,34,47,[50][51][52][53][54]. As concerns computer simulation studies of solutions of simple salts in methanol, we are aware of the following reports [55][56][57][58][59].
Recently, in [60], we explored the properties of NaCl solutions with water-methanol mixed solvent using the SPC/E water [38], united atom methanol [51] and NaCl force field due to Dang [61]. However, the limit of methanolic NaCl solution has not been explored profoundly. Therefore, in the present work we would like to study this particular case more in detail. In addition, our intention is to provide insights into the dependence of the results on the choice of models for methanol and for ionic solutes. This issue has not been investigated so far, in contrast to pure methanol [33,34]. In essence, this is the primary objective of the present contribution. The isobaric-isothermal molecular dynamics computer simulations represent our tools.

Models and simulation details
In the present work, we apply four united atom type models for methanol described in table 1. The CH 3 group is treated as a single site denominated as C. The oxygen and hydroxyl hydrogen sites are denoted as O and H, respectively. More specifically, the methanol model from [52], the model developed in the laboratory of Vrabec [53], the recently developed model from the laboratory of Vega [33], and the model from the TraPPE database [54] are used. They are denominated as L1, L2, OPLS/2016, according to table 1 of [33], and as TraPPE, respectively. The bond lengths and the bond angle for each model are given in table 1.
In table 2 we present the force fields for NaCl employed in the present work. Namely, we use the Joung-Cheatham (JC) [35] and Smith and Dang (SD) [62] models, in close similarity to the description of NaCl aqueous solutions [43]. The values of parameters for ions are taken from JC/SPC/E and SD/SPC/E design, see supplementary material for [43], whereas the crossed, ion-methanol, interaction parameters are obtained using the Lorentz-Berthelot (LB) combination rules. In essence, the three-site SPC/E water model [38] in our procedure is substituted by the methanol models with three sites. We are aware that  the OPLS-2016 force field is designed to operate with the geometric combination rules (CR3 according to the GROMACS nomenclature). In this work, however, we applied solely the LB rules (CR2) in order to have a single and common benchmark for all the models.
In addition, we complement the L2 methanol model [53] with the model for NaCl originally proposed in the work from the same laboratory [63] and later adjusted and applied to describe the properties of aqueous and methanolic NaCl solutions in [57,64]. Again, the LB combination rules are used for this model. In all these models, the ion valency is assumed equal to unity.
Our calculations were performed in the isothermal-isobaric (N PT) ensemble at 1 bar, and at a temperature of 298.15 K. We used the GROMACS software package [65], version 5.1.2.
As concerns the treatment of interactions for each model in question, we would like to mention the following issues. Namely, the non-bonded 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 the precision 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 methanol molecules, the LINCS algorithm was used.
On the other hand, concerning the procedure, a periodic cubic simulation box was set up for each system. The GROMACS 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 energy minimization using the steepest descent algorithm implemented in the GROMACS package. Minimization was followed by a 50 ps N PT equilibration run at 298.15 K and 1 bar using the timestep 0.25 fs. We applied the Berendsen thermostat and barostat with τ T = 1 ps and τ P = 1 ps during equilibration. Constant value of 1.2×10 −4 bar −1 for the compressibility of the solutions was employed. The V-rescale thermostat and Parrinello-Rahman barostat with τ T = 0.5 ps and τ P = 2.0 ps and the time step 2 fs were used during production runs.
Statistics for various properties and various ion concentrations was collected from a set of several N PT runs with duration of 10-20 ns, each started from the last configuration of the preceding run. The time extension for each series of calculations will be mentioned below in appropriate places, if necessary. The number of methanol solvent molecules was fixed at 3000, the number of NaCl molecules (pairs of ions, denoted as N ip ) was in the interval from 0 to 30.

Density of solutions
Our first series of simulations are focused in the description of the dependence of density of NaCl methanolic solution on solute molality, m, in mol/kg units. The simulation results for the entire set or models of this study are compared with experimental data from [66], in figure 1 (a). Slightly less comprehensive, experimental results concerning the density of the system in question at room temperature can be found also in [57,[67][68][69]. The experimental solubility limit has been established at m = 0.238, [70,71].
Inspection of the results shown in figure 1 (a) leads to the following conclusions. An overall location of the lines resulting from different models is determined by the quality of description of the density of pure methanol. All the models correctly describe trends of augmenting solution density on molality. The growth seems to be almost linear. The SD-L1, JC-L1, VR-L2 models overestimate the density of solution at different molalities whereas the JC-TraPPE model underestimates the density. The best agreement with the experimental data is provided by the JC-OPLS/2016 model. A comparison of simulation results for JC-L1 results with SD-L1 data illustrates the role of ions force field. The effect of changing ions force field is small. Evidently, at low concentration of salt, the interactions between ions play minor role. Their effect is slightly more pronounced at a high molality. Additional insight into the quality of presented results can be obtained from figure 1 (b), in terms of reduced density on molality, i.e. the density of solution normalized by the methanol density (this type of plot has been used in [57], cf. figure 5 of that work). In figure 1 (b), we show the results for two models only that exhibit largest difference, the curves from other models flow together, or say lay between the two curves shown in the figure. All the models in question perform well, but the JC-OPLS/2016 predictions seem to be the best. Note that the scale of figure 1 (b) is finer and more adequate, in comparison to the panel describing NaCl solution in figure 5 of [57]. The points of the curves from simulations at a highest molality apparently should correspond to the over-saturated solutions according to the experiment data (theoretical solubility of this salt in methanol has not been established so far up to our best knowledge). However, analysis of the particles configurations by using the vmd software at these conditions (N ip = 30) does not show peculiar features just demonstrating uniformly distributed ions throughout the simulation box.
Therefore, it is of interest to complement these insights by exploring changes of the microscopic structure of NaCl methanolic solutions upon increasing molality for different models in question.

Pair distribution functions, coordination numbers and hydrogen bonding
Evolution of the microscopic structure in the present work is described in terms of various pair distribution functions (PDFs). The pair distribution functions of methanol species have been discussed for various occasions. In order to avoid unnecessary repetition, we would like to comment solely a few specific issues. A set of PDFs describing the structure of pure methanol for three models, the OPLS/2016, TraPPE and L2, is given in figure 2. The curves for each kind of functions resulting from different models are similar, the differences being observed principally in the values of the first maxima, location of the first minima and the phase of the following decaying oscillations. However, these, apparently, subtle differences are manifested in more pronounced discrepancies for various thermodynamic and dynamic properties.
In order to evaluate the quality of predictions of the microscopic structure we resort to the simulated PDFs for the L2 model in [53] (see figure 10). The PDFs from [53] have been compared with their counterparts extracted from the neutron diffraction experiments and empirical potential structure refinement procedure by Yamaguchi et al. [72,73]. The quality of the agreement appeared to be very good. Therefore, the PDFs for L2 model can be considered as a useful reference. The results from our simulations of the L2 model are given by blue lines and symbols in four panels of figure 2. The curves coming from the TraPPE model are very close to the results for L2 model. Small differences can be observed only concerning the height of the first maximum for O-O, O-H, H-H and C-C distributions. On the other hand, the functions coming from the OPLS/2016 model essentially differ in predictions of the values of the first maxima and minima and the phase of further oscillation compared to L2 and TraPPE. It is worth mentioning, however, that the OPLS/2016 predictions are close to the results coming from PHH methanol model [47], as it follows from the inspection of figure 5 shown in [31].
The PDFs describing the distribution of solvent species around cation and anion, from different models of the present study, is given in figure 3 for the lowest molality studied, m = 0.0208. For these functions, the differences between the predictions from different models are more pronounced compared to the previous figure. describe mutual coordination of the species, we use the notion of the running coordination number as common, where ρ j is the number density of species j. The first coordination number, for example, is obtained by putting R = r min , i.e., at the first minimum of the corresponding pair distribution function. As it follows from the g Na-O (r) PDF, the first coordination number for Na + is: 5.5 (VR-L2), 5.7 (JC-TraPPE) and 5.9 (JC-OPLS/2016), respectively. The number of methanol solvent molecules in the second shell is obtained by resting the first coordination number from the second coordination number. It is as follows: 5.4 (VR-L2), 5.6 (JC-TraPPE) and 6.5 (JC-OPLS/2016), respectively. Our result for the VR-L2 model for the first coordination number compares well with the value 5.3 reported for Na + recently at molality m ≈ 0.0313 (cf. table 9 of [57]). On the other hand, the solvation numbers 5.8 and 6.2, for the first coordination number and the number of methanols in the second coordination shell, have been reported in [55] for the ab initio type NaCl solute in PHH methanol for 0.6 molal solution. Moreover, the first coordination number reported here for three models is in agreement with the estimate of 5.7 deduced from IR overtone spectroscopy data in [74] and discussed previously in [55].
The solvation of Cl − anion in methanol can be interpreted by using the Cl-O, Cl-H and Cl-C distribution functions from figure 3. Only the first coordination shell of Cl − is well pronounced, in contrast to Na + . Again, the JC-TraPPE and JC-OPLS/2016 models behave quite similar, just yielding slightly different heights of the first maxima and a different phase of decaying oscillations. The VR-L2 model apparently predicts a more distant position of methanol molecules around the anion. In quantitative terms, these features lead to the following first coordination numbers from the g Cl-O (r) PDF: 5.5 (JC-TraPPE), 5.9 (JC-OPLS/2016) and 5.7 (VR-L2). These estimates are in agreement with VR-L2 data yielding 5.5 (cf. table 9 of [57]). By contrast, the Cl − solvation with PHH methanol model yields a higher coordination number of 7.2, [55], compared to the models of the present study and to the spectroscopy estimates ca. 4, from [74]. From the analysis of ion-methanol atoms distances according to the positions of the first maxima of various PDFs in conjunction, one can construct a schematic picture of cation and anion solvation by methanol. Two examples of this kind of plots are given in figure 4 of [55] and in figure 1 of [58]. The models described here in figure 3 provide a similar picture differing only in subtle details.
Solvation of ions by methanol molecules perturbs the solvent structure and leads to changes of the network of the hydrogen bonds. Hydrogen bonding in pure methanol was explored on several occasions, see e.g., [75] figure 4 (b). This trend preserves in the entire interval of molality studied in the present work. As expected, the n HB decreases with an increasing molality, because the solvation of ions restricts the configurations that satisfy the bonding criteria. Alternatively to the geometric interpretation, the hydrogen bonding can be explored in fine details using energetic criteria, see e.g., [76].
Our final remarks w.r.t. the behaviour of the pair distribution functions concern the distribution of solute ions in the system at a high molality, m = 0.2081. Actually, this value is close to the experimental  It is important to mention that the maxima of g Na-Cl (r) for all models in question are very much higher, in comparison to the aqueous solution of this salt, even close to saturation, cf. figure 4 of [44] for the recently designed so-called Madrid model. Principally, this is the result of strong electrostatic attraction between ions in methanolic solutions due to a low dielectric constant of methanol compared to water. Still, the fraction of CIP is much smaller than that of SSIP, as it follows from the coordination numbers, indicating competition of solvation of ions with the direct Coulomb attraction. The pair distribution functions of the like ions, see figure 5 (b), slowly tend to unity reaching it beyond r ≈ 0.75 for VR-L2 and JC-TraPPE models. The growth of g Na-Na (r) is even much slower within the JC-OPLS/2016 model. This behaviour results from the assumed value of molality and from strong electrostatic repulsion. There is no evidence for anti-phase oscillations or alternating charge distribution for ions and, consequently, for the formation of ion clusters involving more ions than pairs in these solutions. Similar observations concerning the absence of anti-phase oscillations were discussed in [44], while describing Madrid model for NaCl aqueous solutions close to saturation limit. To summarize the above discussion, we would like to note that the models illustrated in terms of pair distribution functions lead to qualitatively correct and physically well grounded conclusions concerning the microscopic structure, in accordance with the previous developments of other authors using computer simulation methods and experimental observations. However, any kind of improvement could be attempted, if the experimental structure factors of such complex systems were available, then fitting of the computer simulation results to, e.g., diffraction experiments data could be performed along the lines proposed in [77,78] for aqueous electrolyte solutions.

Self-diffusion coefficients of methanol and ions
Evolution of the microscopic structure of solutions with molality results in the changes of the dynamic properties. We restrict our attention solely to the evaluation of the self-diffusion coefficients of the species. If the number of particles is sufficiently big to provide good statistics, it is reasonable to obtain the self-diffusion coefficient from the mean-square displacement (MSD) of a particle via Einstein relation, where α is the species index, and the average is taken over all particles (indexed by i) and time origins τ.
To begin with, we used this procedure to obtain, D met . Default settings of GROMACS were used for the separation of the time origins. Moreover, the fitting interval was chosen to be from 10% to less than 50% of the entire trajectory (with time extension not less than 100 ns). The finite size correction has not been taken into account, however. The dependence of MSD on time for systems at different molality, m, for a single model (JC-OPLS/2016) is shown in figure 6 (a). Each MSD line in figure 6 (a) looks perfectly linear, such that the error margin at different m, is small. The entire set of data for different models is given in figure 6 (b). The initial points of curves in figure 6 (b) correspond to pure methanol. Our data agree with the previously reported results, see e.g., table VII in [33] and [22,79]. The experimental result for pure methanol is at 2.44 × 10 −5 cm 2 /s. The curves definitely predict a slowing down of methanol molecules in solutions with an increasing molality. It is difficult to establish how good is the inclination of each curve without experimental data that are unavailable, up to our best knowledge. On the other hand, we are unaware of a similar set of data from molecular dynamics simulations, in contrast to, for example, figure 13 of [44] for aqueous solutions of NaCl.
In order to get more confidence into the quality of the results obtained, we have attempted to perform calculations of D met by an alternative method. Namely, the translational self-diffusion coefficients can be determined from the velocity auto-correlation functions (VACFs) using the Green-Kubo equation, where C vv α (t) is the normalized velocity autocorrelation function of species α, where the average is taken over all particles indexed by i and time origins τ, similar to Einstein formula. To realize the entire procedure, first, we have used the final configuration file from the N PT simulations and performed an additional NVT run for each model with the duration of 100 ps, the time step was chosen at 1 fs. The normalized VACF for methanol center of mass, following from the GROMACS utility, are plotted in figure 7 (a). It appears that the results for L1 and TraPPE models are very similar. On the other hand, the results for L2 just differ from the ones for L1 and TraPPE models by the shape of a shoulder. Finally, the VACF describing the OPLS/2016 model differ from the other curves by the shape of a shoulder, by the depth and position of the principal minimum. An overall shape of the VACF reported here favourably agrees with the result given in figure 4 of [58]. On the other hand, the PHH model for methanol is characterized by a slightly less prominent shoulder, cf. figure 13 of [55]. The Fourier transform of the VACF yields the vibrational power spectrum plotted in figure 7 (b). Two relevant bands of the spectrum were discussed for methanol previously [79]. Namely, the band at frequency around 50 cm −1 and a shoulder at frequency around 150 cm −1 . The first maximum can be assigned to the motion of the particles inside the cage formed by their neighbors whereas the shoulder 23602-10 is attributed to the presence of hydrogen bonds between methanol molecules. From the inspection of the spectrum calculated for different models we see that the shoulder is less pronounced for the OPLS/2016 model, which is due to a lower average number of H-bonds per molecule, in comparison with the other models of methanol, cf. figure 4 (b).
The knowledge of the VACF permits to obtain self-diffusion coefficients of the species. The integral entering Green-Kubo equation should be evaluated as precisely as possible, because it is further multiplied by a big number, k B T/m α for each species. In the case of methanol molecules, we have a good guide to do that using the MSD results. In figure 7 (c), we plot the horizontal lines resulting from the MSD data for each model that should correspond to the integrated normalized VACF. The curves in that figure result from the VACF integration. One needs to find a plateau region, estimate the integral and obtain D met . The plateau is located at different time intervals for different models, as intuitively expected. This kind of procedure has been realized for a set of desired molality points for each model. The results are given in figure 6 (b) as stars joined by dotted lines. In the case of pure methanol, two methods produce quite similar results and the margin of error is less than 2%. On the other hand, trends of slowing down selfdiffusion of methanol molecules are reproduced as well. The best agreement between the data resulting from two methods has been obtained for JC-L1 and JC-TraPPE models. It is also quite satisfactory for the JC-OPLS/2016 model. The inaccuracy of data slightly increases with an increasing ion concentration, the margin of error becomes of the order of 3-4% at the highest concentration studied. The MSD results are apparently more confiable.
Concerning the evaluation of the self-diffusion coefficients of Na + and Cl − ions in methanolic solutions, we would like to make the following comments. We are aware of only scarce data for this specific system. By contrast, there were several reports that describe the dependence of self-diffusion coefficients of ions on molality for aqueous solutions, see e.g., [44,[80][81][82]. Both methods, i.e., the Einstein relation and Green-Kubo equation, were used in above cited works.
Recall, that the number of ion pairs, N ip , in the studied systems is quite low. Therefore, using the mean square displacement procedure on a long time scale leads to an augmenting uncertainty due to sampling statistics, as discussed for example in [80]. In close similarly to that work, we took the slope of the MSD in the interval between 100 and 300 ps to obtain rough estimates for D Cl and D Na through the Einstein relation. The results are collected in table 3. On the other hand, the Green-Kubo equation was used. The VACFs for ions for two models in question in solutions at a rather low molality are shown in figure 8 for illustrative purposes. The curves agree well with the ones reported in figure 2 of [58]. The PHH model for methanol is characterized by a slightly different shape of the VACF for Na + , cf. figure [55]. The values for D Cl and D Na resulting from the VACFs route at a low and at a high molality are given in table 3 as well. Apparently, the agreement between the data coming from two procedures is less pronounced for the JC-OPLS/2016 model. General trend of behaviour of D Cl and D Na is that they decrease in magnitude with augmenting molality in qualitative similarity to what was observed for aqueous solutions, see e.g., figure 14 of [44]. Our data in table 3 for VR-L2 differ but do not contradict the results for NaCl salt in methanol from the report of [57] at molality m ≈ 0.0313. The self-diffusion coefficient for Na + is lower than for Cl − according to the predictions of all models studied.

Solvent dielectric constant
In this subsection we would like to comment on a few peculiar issues of the dielectric response of NaCl solutions with methanol solvent. Since an electrolyte solution is a conductor, the mobile ions move in response to an imposed static electric field. 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 [83]. In several works that describe the application of molecular dynamics computer simulations to electrolyte solutions, the counterpart of the dielectric constant is estimated in terms of the time-average of the fluctuations of the total dipole moment of the system [84,85], obtained by the summation of the dipole moment vectors of the entire set of solvent molecules in the simulation box, where V is the simulation cell volume.
As discussed in detail in [86] for aqueous electrolyte solutions with completely dissociated solute molecules (NaCl in water is an example, see table 7.6 of [86]), the static permittivity of the sample, ε, decreases with increasing the electrolyte concentration. This permittivity represents the static permittivity of the solvent in the solution. Moreover, the concentration dependent decrease of the solvent permittivity defines the dielectric decrement. For aqueous solutions of NaCl it was discussed using, for example, the micro-field approach [87] and the results are in favourable agreement with experimental data [88]. We would like to briefly report this property as a function of ions concentration for different models of solutions. Namely, in figures 9 (a) and 9 (b), the time evolution of ε solvent at different concentrations of ions in a solution, is shown for the JC-OPLS/2016 and VR-L2 models, respectively. It is known already that quite long MD runs are necessary to adequately describe the static dielectric constant of pure insulating liquids. Even after 100-120 ns runs, the data in these two panels do not show convergence to a perfectly well defined constant. Therefore, the averaging over the final time interval was performed. However, it is important to mention that the static dielectric constant of methanol predicted by the OPLS/2016 model is closer to the experimental value equal to 33.1, compared to all other models of this study. In addition, one can observe that the solvent dielectric constant decreases with an increasing ion concentration, figure 9 (c). It is impossible to evaluate how good is the magnitude of the decrement of curves without experimental data for this particular property. As it follows from the data in figure 9 (c), the downward inclination of the curves is less pronounced, in comparison with the recent report concerning the modelling of the permittivity of electrolyte solutions [89] (see figure 5 of this work for the dependence of ε solvent on molarity for various salts in methanol).

Surface tension
The simulations at surface tension calculations aiming at each point of composition axis were performed by taking the final configuration of particles from the isobaric run. Next, the box edge along z-axis was extended by a factor of 3, generating an elongated box with a liquid slab and two liquidmixture-vacuum interfaces in the x−y plane, in close similarity to the procedure used in [90]. The total number of methanol molecules, 3 × 10 3 , is reasonable to yield an area of the x−y face of the liquid slab sufficiently big. The elongation of the liquid slab along z-axis is satisfactory as well. The executable molecular dynamics file was modified by deleting a fixed pressure condition just preserving V-rescale thermostatting with the same parameters as in N PT runs. Other corrections were not employed. The values for the surface tension, γ, follow from the combination of the time averages for the components of the pressure tensor, where P i j are the components of the pressure tensor along i, j axes, and . . . of 10 ns, and obtained the result for γ by taking the block average. However, only the results coming from the blocks yielding most close P xx and P y y were taken into account. The results are given in figure 10. Our simulation data for pure methanol are shifted downward or say, systematically lower than the ones reported in [33], presumably because of the difference in technical details of simulations, e.g., the geometric combination rule (CR3 in GROMACS nomenclature) was used for OPLS/2016 model in [33]. The experimental result for pure methanol is at γ = 22.51 mN/m. However, the sequence of the magnitude of the results for different models in the present work and in [33] is the same. The surface tension values definitely increase with an increasing ion concentration. This trend is not strong, however. From the density profiles of methanol species (the OPLS/2016 model is chosen just for illustration) at m = 0 and at m = 0.3121, we conclude that the interface becomes slightly sharper upon adding ions to methanol [ figure 11 (a)]. The density of the liquid slab is lower at both molalities than for the bulk system, cf. figure 1 (a). Apparently, if the number of particles in simulations increased, the interface width could slightly decrease and the surface tension values might become higher. On the other hand, as it follows from figure 11 (b), the interface layer is constituted by methanol molecules. There is no preferential adsorption of ions at the interface, as expected for the nonpolarizable methanol and ions models from the macroscopic electrostatic arguments. The ions are expelled from the interface as witnessed by the density profiles of ions in figure 11 (b). This issue for aqueous solutions of salts was discussed in detail in [91].

Summary and conclusions
We have performed a rather extensive set of molecular dynamic simulations to study various properties of NaCl solutions in methanol at room temperature and ambient pressure. The models for solutions involve united atom type description for methanol using L1, L2, TraPPE and OPLS/2016 force fields and the JC, SD and VR models for ions. The combination of models of the solvent and solute species has been inspired by the previous developments for aqueous solutions of NaCl.
We explored the evolution of density, surface tension, solvent dielectric constant and self-diffusion coefficients of the species with an augmenting amount of NaCl solute in the solution. The microscopic structure is described in terms of the pair distribution functions, the first coordination numbers of species and the number of the ion pairs. Trends of behaviour of different properties of the system are determined by the balance of effects resulting from non-electrostatic and electrostatic interactions. This is manifested, for example, in a very different solubility of this salt in methanol solvent compared to water. From a comparison with a limited set of experimental data for the systems in question and with the results of other authors on related systems, we can conclude that the predictions obtained are qualitatively correct and give a physically sound picture of the properties explored. It seems that the OPLS/2016 model for methanol combined with JC force field for ions performs best, compared to other combinations of the force fields.
At the present stage of the development, the missing and very important elements worth a much more detailed investigation are numerous. 1) Within the employed modelling, it is worth to explore the shear viscosity and electric conductivity, as it has been done for NaCl solutions in e.g., [30,43,44,57]. 2) It is worth to explore the extension of the present models by involving either four-site or all-atom description of methanol combined with different force fields for ions. 3) Ampler insights into the trends of behaviour of various properties of this kind of solutions would follow from the study of a set of alkali halide salts in the same solvent. 4) The role of combination rules should be explored as well, see e.g., [92] for the case of aqueous solutions of NaCl.
The most challenging and remaining task, however, is to establish the solubility limit for the models in question by using the tools tested for NaCl aqueous solutions in [43,44,93]. We hope to address some of these issues in a future work.