Towards the correct microscopic structure of aqueous CsCl solutions with a comparison of classical interatomic potential models

The structure of aqueous CsCl solutions was investigated by classical molecular dynamics simulations (MD) at three salt concentrations (1.5, 7.5, and 15 mol %). Thirty interatomic potential sets, based on the 12-6 Lennard-Jones model, parametrized for non-polarizable water solvent molecules were collected and tested. Some basic properties, such as density, static dielectric constant, and self-diffusion coefficients, predicted by the force fields (FF), were compared with available experimental data. The simulated particle configurations were used to calculate the partial radial distribution functions (PRDF) and the neutron and X-ray total scattering structure factors (TSSF). The TSSFs were compared with experimental data from the literature, to find the best FF models, which describe the structure correctly. It was found that, though several of the thirty models failed in the tests, some models are compatible with the measured data. Values of the structural parameters consistent with the experiments were determined (such as water-ion distances, the average number of water molecules around the ions, average number, and distance between anion-cation contact ion pairs, water-water hydrogen bonds). It was shown, that in addition to models in which the number of contact ion pairs is too high, models in which this number is too low are also unable to reproduce the experimental data.


Introduction
Simulation techniques play an essential role in studying the properties of molecular liquids. The same applies to describing the structure of aqueous salt solutions. Since at least four constituting elements (O and H from water, anion and cation from the salt) are present in such solutions, there are ten partial radial distribution functions (PRDF). There should be ten independent (diffraction) experimental data sets to determine all the ten PRDFs, which is impossible. Thus the use of simulation techniques is necessary to determine the structure.
Computer simulation techniques have been applied for investigating molecular liquids for nearly 70 years [1]. Initially, Metropolis Monte Carlo (MC), classical molecular dynamics (MD), and reverse Monte Carlo (RMC) methods yielded the most publications on this topic; but today, MD simulations using more complex potential models and ab initio molecular dynamics (AIMD) simulations have come to the fore. However, the classic MD simulations, using rigid, non-polarizable models, are still widely used nowadays due to their widespread acquaintance and lower computational costs. Moreover, the development of such interatomic potential parameters also continues (see e. g., Ref. [2]).
The success of the classical MD simulations depends fundamentally on the interatomic potential model (force field, FF) chosen. It was shown for aqueous NaCl solutions that there are crucial differences in the performance of the force fields, even in similar types of models [3]. Some of the tested potential models gave false results even for basic properties such as the solubility range: precipitation was observed in the simulated system at concentrations below the (experimental) solubility limit. It was found in a recent paper that pairwise additive MD models give poor agreement with experimental EXAFS data on aqueous CsCl solutions [4]. The authors of that paper concluded that it is necessary to include many-body effects in the simulations to get the correct description of the structure. Since the usability of the classical models highly depends on the FF parameters, the validity of such a statement should be tested against more FFs, rather than just the one used in Ref. [4].
Our previous study investigated the applicability of 29 different 12-6 type Lennard-Jones potential models for highly concentrated aqueous solutions of LiCl [5]. The comparison of several FFs with experimental diffraction data made it possible to describe the structure of aqueous LiCl solutions: the hydrate sphere of the ions (number and distance of water molecules around the ions), as well as the presence of contact and solvent separated ion pairs [6]. This technique has been completed with the analysis of the hydrogen-bonded (H-bond) network of water molecules and chloride ions [7], showing that considering the halide anions to be part of the H-bonded network makes a more meaningful interpretation of the structural details and is also helpful to separate better and worse performing potential models.
The structure of CsCl solutions and the hydrate sphere of Cs + and Clions have been studied since 1967 [8] applying experimental and simulation techniques. Reviews of previous results are collected in Refs.
[9 -11]. Experimental techniques were applied in several papers, e. g.: X-ray diffraction (XRD) [ Some of the papers combined experimental data and simulation techniques, such as in Refs. [4,12,13,17], which helps to obtain more realistic results.
According to the papers mentioned above, the hydrate sphere of Cs + ions consists of 6 -10 water molecules, with a Cs + -O distance of about 3 -3.3 Å. Around the Clanions, 5 -7 water molecules can be found with Cl --O distances of about 3.1 -3.2 Å. The number of hydrating water molecules decreases while the formation of contact ion pairs increases with increasing salt concentration. The Cs + -Clcoordination number is around 2 -3 in the most concentrated solutions, with a Cs + -Cldistance around 3.35 -3.5 Å.
The values of the structural parameters found in the literature have a fairly broad distribution, either obtained by simulation or experimentally. Since the structural parameters cannot be measured directly, the interpretation of the measured data also influences the (so-called) experimental values. This also increases the uncertainty of these values. During the development of FFs, the appropriateness of targeted experimental data is crucial; hence it is necessary to provide more accurate target values.
Most FFs consist of complete parameter sets for all halogen and alkaline ions. The target properties of these FFs are different, but they are usually not the structural parameters of concentrated solutions. The ability of one model to describe the structure of one type of salt solution does not mean that the same FF will give correct values in the case of another solution. (It will be shown later, that FFs, which predict the structure of aqueous LiCl solution well, can completely fail for aqueous CsCl solution and vice versa).
30 different 12 -6 Lennard-Jones type ionic interatomic potential models, developed for rigid, nonpolarizable water models, are investigated in this study. They are used to simulate highly concentrated aqueous CsCl solutions. The appropriateness of the models is determined by comparing the simulated results with experimental data, concerning both the structure and some basic properties. PRDFs and total scattering structure factors (TSSF) are calculated from the obtained atomic configurations and are compared with diffraction data. The structural results of the best models are further investigated to describe the structure in more detail.
The purpose of this study is twofold: (1) to ascertain which FF parameters are appropriate for the description of aqueous CsCl solutions and (2) to determine the structural parameters (bond lengths, coordination numbers) more precisely to get values, which can be targets for future FF-development.

Investigated concentrations, experimental ND and XRD diffraction data
Aqueous cesium chloride solutions were investigated at three concentrations, namely at 1.5, 7.5, and 15 mol %. Neutron and X-ray diffraction data are known for all these concentrations [12]. The exact concentrations (number of ions per water molecules) and the corresponding number densities were taken from Ref [12] and are shown in Table 1. Simulations were performed in cubic simulation boxes, using periodic boundary conditions, the box sizes were calculated from the number densities. The number of atoms in the boxes was around 10000; the exact number of ions and water molecules are summarized in Table 1.

Interatomic potentials (force field models)
Pairwise additive, non-polarizable potential models were investigated, in which the intermolecular interactions between atoms i and j are described by the 12-6 Lennard-Jones (LJ) interaction and the Coulomb potential: (1) Here rij is the distance between two particles, i and j, qi and qj are the point charges of the two particles, and ε0 is the vacuum permittivity. εij and σij represent the energy and distance parameters of the LJ potential. The ionic εii, σii, and qi values for the investigated models are collected in Table 2 (also shown in Fig. S1 in the Supplementary Material (SM)). The tested FFs apply one of the following water models: SPC/E [28], TIP3P [29], TIP4P [29], TIP4P-Ew [30] and TIP4P/2005 [31] (also shown in Table 2). All of these water models are rigid and non-polarizable, the first two are 3-site models (3 point-like atoms representing the oxygen and the two hydrogen atoms), the other three models use a fourth, virtual site also (the partial charge of the oxygen atom is located on the virtual site). The LJ parameters of the hydrogen atom are 0 for all of these water models; the other parameters are collected in Table 3.
The LJ parameters between unlike atoms (anion-cation and ion-oxygen) are calculated according to combination rules. Most tested models (except for one) use one of the following rules: a) the geometric combination rule (where both εij and σij are calculated as the geometric average of the homoatomic parameters) or b) the Lorentz-Berthelot combination rule (in which εij is calculated as geometric, while σij as the arithmetic average of the relevant parameters). The applied combination rules are shown in Table 2.
The tested interatomic potential sets are listed below. Several of the FFs studied here were also tested for aqueous LiCl solutions; a more detailed description of those FFs can be found in Ref. [5].

Ion parameters with SPC/E water:
 Dang-Smith parameter set (DS) [32,33]: One of the most frequently used FFs, it was applied in Ref.
 Joung-Cheatham III set (JC-S) [34]: Parameters were fitted for hydration free energies, lattice energies, and the lattice constants of salt crystals. The authors developed parameters for three different water models (all of them are tested here).
The ion-oxygen parameters were determined, leaving the choice of the combination rule free.
Therefore geometric and LB rules were also tested here.  [39,40]: These parameters were determined to reproduce the experimental free energy of solvation for a given single ion, the ion-pair properties were taken into account through fitting osmotic coefficients. For the cesium ion two different LJparameters are given, which were similarly appropriate in the criteria mentioned above; the notation of Ref. [39] is followed here, the numbers 6 and 9 are used to distinguish the two sets.
of charge-scaling was tested with two previously presented models: RL (with SPC/E water) and Li-IOD-T (with TIP4P-Ew water) models.

Molecular dynamics simulations
Classical molecular dynamics simulations were performed by the GROMACS software (version 2020.6) [49]. Ions and water molecules were initially placed into the boxes randomly. Water molecules were kept together rigidly by the SETTLE algorithm [50]. Van der Waals interactions were truncated at 10 Å, with added long-range corrections to energy and pressure [51]. Coulomb interactions were treated by the smoothed particle-mesh Ewald (SPME) method, using a 10 Å cutoff in the direct space [52,53]. Energy minimization was carried out using the steepest descent method. After that, the leapfrog algorithm was used for integrating Newton's equations of motion, using a 2 fs time step.
Equilibration was carried out at constant volume and temperature (NVT) at T = 300 K. In the first equilibration stage, which was 200 ps long, the Berendsen thermostat [54] was used with a time constant τT = 0.1, after that the Nose-Hoover thermostat [55,56] was applied to keep the temperature constant, with τT = 2.0. After a 4 ns long second equilibration stage, configurations were saved in every 2 ps during a 20 ns long production run. For density determination, a 4 ns long NpT simulation was carried out. The pressure was kept at p = 10 5 Pa by Parrinello-Rahman barostat [57,58], using a τp = 2.0 coupling constant. In these simulations the Nose-Hoover thermostat was used with τT = 0.5 coupling.

Calculation of the partial radial distribution functions and the structure factors
Every 10 th (altogether 1000) configurations were used for calculating the partial radial distribution functions (PRDF, gij (r)) with the 'gmx rdf' program of the GROMACS package. The partial structure factors (Sij (Q)) were calculated as Here Q is the amplitude of the scattering vector, and ρ0 is the average number density. The neutron and X-ray total scattering structure factors (TSSF) were obtained as (3) Here wij X,N are the scattering weights for X-ray (wij X ) and neutron (wij N ), they are given by equations (4 and 5): where δij is the Kronecker delta, ci denotes the atomic concentration, bi is the coherent neutron Here Smod and Sexp are the model and experimental structure factors, and the summation is over the Qi experimental points.

Hydrogen bond calculations
Every 100 th (altogether 101) configurations were used for hydrogen bond analysis. Hydrogen bonds (HBs) were identified using energetic criteria: two water molecules are identified as H-bonded if the intermolecular distance between an oxygen and a hydrogen atom is less than 2.5 Å, and the interaction energy is less than -12 kJ/mol. A similar definition is used in the case of the chloride ion using the same energy boundary and a cut-off distance (H...Cl -), which is defined as the first minimum value of the Another possibility offered by the detailed comparison of a multitude of FFs is to analyze subtle variations between members of model families in order to identify trends. Beside the effect of charge scaling, which was detailed previously, such example is the family of the Horinek models (HS, HM and HL). The analysis of these models show that the density and water self-diffusion values does not vary significantly between the HL-HM-HS models, but the static dielectric constant increases as the interaction strength decreases. At the same time the ionic mobility decreases and the tendency to precipitate increases.

Neutron and X-ray total scattering structure factors
Partial radial distribution functions were calculated from the obtained trajectories, they will be discussed in detail in the next section. Neutron and X-ray TSSFs were calculated from the PRDFs, and are compared with experimental results from Ref. In the case of the less concentrated composition (1.5 mol %), S N (Q) depends mainly on the applied water model: FFs using 4-site water models perform best, while the agreement with experimental data is the poorest for models using the TIP3P water. The R-factors for the FFs with 4-site water models are around 10 %, they are 16-18 % for the FFs with the SPC/E model, and more than 20 % for the models with the TIP3P water. At higher concentrations, the applied ionic parameters have more influence on the ND fits, and the R-factor has broader distribution: its value is between 10 % and 50 % for the 15 mol % sample. In general, the FFs with 4-site water models give better results. However, there are also significant differences between these FFs (see, for example, Li-HFE-T and Li-IOD-T models).
Concerning the FFs with 3-site water models, the role of the water model seems to be less important: there are huge differences between the FFs: see e. g. RH (R = 52 %) and RL (R = 17 %) models. The most significant deviations are at around Q = 1 Å -1 and Q < 1 Å -1 . In the experimental curve, a small pre-peak emerges around Q = 1 Å -1 , while in the simulated curves, this peak is missing for some FFs The R-factor values of the X-ray TSSFs (S X (Q) functions) obtained by simulations using different FFs are between 6% and 20% in the case of the 1.5 mol % solution, except the HS-LB model, which has an R-factor as high as 45 %. In the 15 mol % sample, the R-factors are between 6% and 38%, except the HS-LB model with 80% and the HS-g model with 67%. Beside the HS-LB model, which is definitely precipitated even at the lowest concentration, the HS-g model also has salt precipitations at 15 mol %.
Some other models show clustering of the ions at the highest concentration, e. g. JC-S or FN6 models, which FFs also produce high R-factors. The performance of the Madrid model is surprisingly poor, despite the absence of ion clusters in the configurations, which suggests that the lack of ion clusters is not sufficient, not even necessary to ensure a correct structural model for aqueous CsCl salt solutions.
In some simulated curves, there is a small pre-peak around Q = 1 Å -1 , similarly to the ND curve, however, in the experimental data, only a shoulder can be seen. In the experimental curve, the 1 st (or main) peak is around Q = 2.1 Å -1 . A 2 nd peak (for the 1.5 mol % solution) or a shoulder (for the 15 mol % solution) is at around Q = 2.6 -2.8 Å -1 , followed by a 3 rd peak at about 4 -4.5 Å -1 . Though these peaks are identifiable on the simulated curves, their positions and amplitudes vary from FF to FF.
The ratio of the amplitudes of the 1 st and 2 nd peaks (at Q = 2.1 Å -1 and 2.8 Å -1 ) is also different.
Concerning the 1.5 mol % sample, there are 12 models with similarly good R-factors. As the concentration increases, the R-factor of some models remains nearly constant (e. g. DS, Li-HFE-T), but for some FFs the R-factor increases (e. g. MP-T, RDVH). In some cases, the more concentrated solution has a better fit (mostly for the models with TIP3P water, see e. g. Li-HFE-T3). The lowest and was also tested for experimental EXAFS data on CsCl solutions [4], is not among the best models concerning the structure of aqueous CsCl solutions. Thus the finding of ref.
[4] that classical, nonpolarizable MD models are not suitable for describing the structure of aqueous CsCl solutions is unfounded since there are classical models performing significantly better than the model they have chosen.

Hydration shell of the anions (Cl -)
The gClH(r) and gClO(r) curves for some FFs are shown in Figs. 8 and 9. They have a well-defined 1 st peak, followed by a relatively deep minimum for all investigated FFs. This peak originates from the water molecules hydrating the anion. The positions of the first peaks (rClH and rClO) do not depend on salt concentration in the investigated concentration range, but they depend strongly on the applied ionic potential model, see Table 4 and Fig. 10. The position of the first minimum is well defined and concentration-independent for the Cl --H pairs. In the case of the Cl --O pairs, the first minimum is not as deep, it is indistinct, and slightly depends on concentration (see Tables S6 and S7).
The coordination numbers are calculated from the PRDFs as: where cj is the concentration of the j th species, and R is the upper limit of the coordination sphere.
Usually, the upper limit of the integration is defined as the first minimum after the first peak. The NClH, NClO coordination numbers (calculated up to the first minima) are also collected in Table S6 and S7, and shown in Fig. 11.
The NClH and NClO coordination numbers are close to each other (within 10 %), although NClO is slightly higher than NClH. This shows that water molecules connect to Clions with only one of their H atoms.
The water molecules are oriented with one of their H atoms towards the anion; the other H atoms of these first shell water molecules produce the 2 nd peak of the gClH(r) curve. This orientation is also verified by the Cl --H-O angle, which is nearly 180°, as shown for some models in For most models 5.7 -7 water molecules can be found on the 1 st hydration shell of Clions at 1.5 mol %, which number decreases to 4.8 -6.5 for the 7.5 mol % solution and 3.7 -5.9 in the 15 mol % The first hydration shell of the cesium ions can be studied through the gCsO(r) curves; for some FFs, see  Tables 4 and S8, and Figs. 10 and 11. For most FFs, the NCsO coordination numbers are between 6.6 -9, 6 -8.5, and 5 -8 for the 1.5, 7.5, and 15 mol % solutions, respectively. Some FFs give lower hydration numbers; these are the models with precipitations (HS-LB, HS-g). In contrast, there are two models with significantly higher NCsO values: the RH and RM models, which also have the longest bond distance; in these models, the 1 st peak is strongly asymmetric with an elongated decrease, the limit of the 1 st hydration shell cannot be reliably determined. For the best model, the KS-Li-IOD-T, NCsO is 8.56, 7.58, and 7.4 for 1.5, 7.5, and 15 mol % solutions, respectively.

Contact ion-pairing
The first peak of the gCsCl(r) curve (around 3.5 Å) originates from the Cs + -Clcontact ion pairs ( Fig.   13). It is a sharp, well-defined peak, followed by (for most models) a rather deep minimum (except Madrid FF, see below). The 2 nd peak is small and broad; it is produced by those Cs + -Cldistances, for which the ions are connected through a common water molecule (Cs + -O-H -Clmotifs). The position of the 1 st peak (rCsCl) depends on the FFs, its value is in the 3.08 -3.80 Å range, see Table 4 or The NCsCl coordination number (the number of contact ion pairs) is usually in the range of 0.16 -0.48, 0.68 -1.65, and 1.76 -3 in the 1.5 mol %, 7.5 mol % and 15 mol % solutions, respectively, see Table   S9 or Fig A similar rule holds for the Cs + ions: NCs=NCsCl+NCsO is more or less constant: the difference is mostly below 5 %. The total coordination number of Cs + ions obtained by different FFs has a broad distribution, however, between 7 and 9.5. For the most appropriate models, it is 8.5 < NCs < 9.5 (for the KS-Li-IOD-T, NCs is 9.23).
The Madrid model is unique in this comparison because it was developed specifically to produce an extra-low number of ion pairs. The shape of the gCsCl(r) curve obtained by this FF is different from that of the other models (see Fig. 13). There is no sharp first peak in the 3 -4 Å region. The first visible peak is around 5 Å, with a shoulder on the left side. The peak around 5 Å originates from the solventseparated ion pairs, while the shoulder represents the contact ion pairs. Since the contributions of contact and solvent-separated ion pairs overlap, the number of contact ion pairs cannot be determined PRDFs. The ion-water nearest neighbors were discussed in the previous sections. The water-water, hydrogen-bonded connections can be seen as the first intermolecular peak in the gOH(r), which is located around 1.8 Å (1.76 -1.84 Å); the position of the peak depends only on the applied water model (Fig. 14). The NOH coordination numbers (which are uniformly determined up to 2.4 Å for all FFs) are collected in Table S10 and shown in Fig. 11. As the salt concentration increases, the number of waterion connections increases and the number of H-bonded water pairs decreases. The water molecules connect with their H atoms to O atoms of other water molecules (H-bond) or to Clions. The NH=NHO + NHCl coordination number is around 0.9 -1 for all investigated FFs; the higher numbers belong to the higher concentrations, which means that nearly all of the H atoms of water molecules participate in Hbond or water -Clconnections.
In pure water, the oxygen atoms of water molecules link to two water molecules through H-bonds. In the salt solution, they lose one (or two) of their H-bonded water pairs to form water -Cs + pairs instead.
Since the upper limit of the hydration shell of Cs + ions and thus the NOCs coordination number has high uncertainty, the NO = NOH + NOCs total coordination of O atoms is also quite uncertain. Its value is more or less around 2 for the 1.5 and 7.5 mol % solutions, but it is mostly higher than 2 in the 15 mol % solution. However, the lack of a well-defined minimum in the gCsO(r) can cause the overestimation of the NCsO coordination number.
The first peak in the gOO(r) PRDF (around 2.75 Å) originates mainly from the H-bonded water molecules ( Fig. 15). However, the first minimum of this peak is not well-defined, and for the FFs applying the TIP3P water model, there are no minima at all. Thus the NOO coordination number was determined uniformly up to 3.2 Å for all FFs (Table S11) which are paired to a common ion. (Water molecules, which are in the hydrate sphere of the same ion).

Structural origin of special features in the total scattering structure factors
The shapes of the total scattering structure factors depend on the applied potential model. The TSSFs are calculated from the PRDFs, so dissimilarities in TSSFs result from differences in PRDFs. In this section, two Q regions are investigated to find the reason for the difference in TSSFs obtained by various FFs. Only the results and the main steps of these calculations are presented here; a detailed description can be found in the SM. Two significant differences can be observed between the potential models in the S N (Q) curves of the 15 mol % solution at low Q values: around Q = 1 Å -1 and below that value. At Q = 1 Å -1 , there is a small peak in the experimental curve, it is usually called "pre-peak", or "first sharp diffraction peak" (FSDP).
The presence of this peak is often associated with the medium range order in the system. The existence and height of this peak in the simulated curves highly depend on the applied potential model (see Fig.   4). Below the pre-peak region, as Q → 0, some of the calculated S N (Q) functions begin to increase, which is in contrast to the behavior of the experimental curve (e. g. for the MS, Li-HFE-S, Li-HFE-T, Li-HFE-T3, DVH, RDVH, JC-S, and FN6 models). This behavior was observed previously (see Ref. [12]), in that study, the DS FF parameters were used in the MD simulation. This result also shows that many structural motifs cause the pre-peak of the S(Q) curve together; the presence of the FSDP cannot be assigned to one special motif in the solution.

H-bonded network
Traditionally, hydrogen bonding (H-bond) means the interaction between water molecules; however, an identical mechanism is responsible for the connection between the chloride ion and the hydrating water molecule. The halide ion can be considered a constituent of the hydrogen-bonded network, as was discussed in a previous paper [7]. Here the H-bonded network structure of aqueous CsCl solutions is investigated for some selected models.
The H-bonded water-water and water-chloride ion pairs were specified via an energetic definition of Hbonding, see section 2.5. The size distribution of H-bonded clusters for the water-anion system and the water subsystem (taking into account only the water-water connections) are shown in Figs. 18 and 19.
Considering the water-anion system, the size of the most extended clusters are in the order of the number of water molecules plus Clions for all concentrations; only a slight decrease can be observed as the salt concentration increases. There are small differences between the FFs: the size of the longest clusters are slightly bigger for models in which the ion clustering is less frequent. The network of water molecules is percolating for the 1.5 and 7.5 mol % solutions; however, the size of the clusters decreases significantly as the salt concentration increases. There is a huge difference between the FFs in this regard: those models, in which there are more ion clusters, contain more water-water connections as well, and the water H-bonded network is percolating even for 15 mol % salt content. However, in models, which has a weak ion-paring tendency, the size of the water clusters is significantly smaller than the total number of water molecules.

Conclusions
Aqueous CsCl solutions were investigated by classical MD simulations. Thirty interaction models, all of them using 12-6 LJ potentials and rigid, non-polarizable water molecules, were examined. The calculations were made at three different concentrations from 1.5 mol % to 15 mol %. Density, static dielectric constant, and self-diffusion predictions of simulations were compared to experimental values from the literature. The partial radial distribution functions were calculated for all FFs, together with neutron and X-ray total scattering structure factors that were also compared with previous experimental results.
The simulated density values are close to the experimental values, the static dielectric constant values are smaller than the experimental ones. The self-diffusion coefficient values are close to the experimental ones at the lowest concentration; however, they are significantly smaller at the highest investigated concentration. Two FFs, namely HS-LB and HS-g, precipitate below the solubility limit.
The structure parameters predicted by the potential models were investigated in detail. There are huge differences between the models; however, comparing the calculated structure factors with those obtained in neutron and X-ray diffraction measurements made it possible to find the best performing models. Special features that appear in the structure factors calculated from different models also Some additional comments about the performance of the models: (i) The charge scaling is successful concerning both the physical properties and structural behavior. The density decreases, the static dielectric constant and the self-diffusion coefficients increase as the charge of the ions decreases.
(ii) The suitability of the models for one type of salt does not automatically mean their good performance for another salt.
(iii) Though the number of contact ion pairs should not be too high to avoid precipitation below the solubility limit, the experimental data clearly shows that some amount of such pairs is present in the CsCl solutions.
(iv) Structure parameters (bond length, coordination numbers) should also be taken into account during the development of force field models. Thus it is very important to determine them (using experimental and simulation techniques together) as accurately as possible.
(v) Despite the serious limitations of simple, non-polarizable, point charge classical potential models, they can be suitable to describe the structure of aqueous salt solutions, but it is necessary to use experimental diffraction data to find (or create) the best possible parameter sets for that purpose. [

Density
A 4 ns long NpT simulation was carried out to determine the density predicted by the tested model. The pressure was kept at p = 10 5 Pa by Parrinello-Rahman barostat [S1,S2], using a τp = 2.0 coupling constant. The Nose-Hoover thermostat was used in these simulations, with τT = 0.5 coupling. Densities were calculated every 2 ps from 2 to 4 ns, and the average of the 1000 values was determined.
Density values obtained by MD simulations are collected in Table S1. In Fig. 1

Static dielectric constant
The static dielectric constant was calculated according to Eq. S1: where M is the total dipole moment, V is the volume of the system, kB is the Boltzmann-constant, ε0 is the vacuum permittivity, and T is the temperature. All saved configurations were used as an input for the calculations performed with the 'gmx dipole' program of the GROMACS software. The average of the ε values in the last 10 ns part of the ε(t) curves (10 ns -20 ns) was determined and taken as the final result.
Some typical ε(t) curves are shown in Fig. S4. Long trajectories are necessary until the curves converge to the final ε values, as was previously demonstrated for pure water in Ref [S6]. A 6 ns long trajectory is enough for pure water, as was found in Ref [S6]. However, in concentrated salt solutions, ε(t) curves are not saturated even after 8 ns [S7]. Thus 20 ns long trajectories were used in this work; and the last 10 -20 ns of the curves were averaged to determine the ε value. The 20 ns was sufficient for ε(t) to saturate at all concentrations tested, except for two FFs: HS-LB and HS-g.
ε values obtained from the simulations with different FFs are shown in Fig. 2 of the main text, and Table S2. The experimental values at the given concentration were calculated using the formula of Ref.   The ratio of the static dielectric constant of the solution to that of pure water was also calculated to eliminate the effect of the water model and to focus on the differences between the FFs of the ions. The

Self-diffusion coefficients
The self-diffusion coefficients were calculated from the mean-square displacement (MSD) using the Einstein-relationship (Eq. S2): Here A refers to either Cs + , Clions, or water molecules. The full trajectories were used by restarting the MSD-calculation every 10 ps. The 'gmx msd' program of the GROMACS package was used for the calculations. The linear part of the MSD-t curve was fitted to determine DA. The MSD-t curves of the water molecules were linear for the full 0 -20 ns region; they were fitted in the 2 -18 ns part.
However              In Fig. 4 and 5 (in the main text) the main features of the total scattering structure factors can be seen, the experimental ones together with those obtained from simulations using some selected potential model. In the S N (Q) curve of the 15 mol % solution a peak can be seen about at 1 Å -1 , which is often called as pre-peak, or first sharp diffraction peak. The presence of this peak is often associated with the medium range order in the system. As it can be seen in Fig. 4 (in the main text), the presence and height of this peak in the simulated curves highly depends on the applied potential model. This can help us to investigate the structural origin of this peak. The amplitude of S(Q) at a given Q value is determined as a weighted sum of the partial Sij(Q) curves, which are calculated from the PRDFs (see Eqs. 2, 3 in the main text) The low Q part of the S N (Q) curve and the neutron weighted partial structure factors (wij*Sij(Q)) of some selected FFs are shown in Fig. S18.  (c) H-H number of H-bonded water molecules, the rClH bond lengths as well as the number of hydrating water molecules can affect the shape of the H-H PRDF and cause the pre-peak of S N (Q).

Neutron and X-ray
Another interesting part of the simulated S N (Q) curves is below the pre-peak region. As Q → 0, some of the calculated S N (Q) functions begin to increase in contrast to the experimental curve (e. g. for the MS, Li-HFE-S, Li-HFE-T, Li-HFE-T3, DVH, RDVH, JC-S, and FN6 models). This behavior was observed previously in Ref. [S3] also, in that study the DS parameters were applied in the MD simulation. To investigate the origin of this discrepancies a comparison, similar to that detailed in the previous paragraph, was performed at Q = 0.5 Å -1 . Figure S21. Values of the weighted partial structure factors at Q = 0.5 Å -1 for the selected force field. Note: the y-scales, although starting from different values, span the same range to make the comparison easier. Here the main differences between the FFs are in the H-H and O-H partial structure factors, whose effect is partly compensated by the Cl --H, Cs + -H and Cl --O partials (Fig. S21). In the r-space, several regions contribute to the difference in the 2 -8 Å regime, see Fig. S22. Let's choose the DS and Gee FFs, which are significantly different at Q = 0.5 Å -1 , but are similar at higher Q values. Comparison of the DS and Gee FFs shows, that the main difference between these two FFs is in the 2 -4 Å region both in H-H and O-H PRDFs. As the Cl --H bond distance is similar in these models, the difference is mostly caused by the different number of H-bonded water molecules: the NOH coordination number is 1.23 for DS and 1.1 for Gee model. This difference seems to be moderate, however the shape of the (b) H-H with the decrease of ion-water pairs (the hydration shells of ions are connected for clustered ions, they create a common hydrate shell), while the water molecules can establish H-bonds between themselves.
This phenomenon leads to the increase in the structure factor at low Q values.