The “good,” the “bad,” and the “hidden” in neutron scattering and molecular dynamics of ionic aqueous solutions

The fundamental properties of ion–ion and ion–water interactions regulate various biological and geochemical processes ranging from protein and enzyme activation to desalination and mineral nucleation.^1,2 Calcium is probably the most studied divalent cation due to its abundance and importance. Biological environments possess small amounts of calcium ions resulting in only dozens of millimoles per liter of water. Still, this amount is enough to stimulate or suppress protein function or drive intracellular processes known as calcium signaling.³ Geochemically relevant processes occur at or near the solubility limit, where the formation of ionic clusters can be a precursor to crystallization of, e.g., calcium carbonate, which is the most common calcium-containing compound on Earth.⁴ 

Molecular dynamics (MD) simulations are widely used to study, explain, and even predict such phenomena. Although the solution thermodynamics strongly depends on the concentration, it is common practice that models optimized at infinite dilution are used to investigate the processes over a wide range of thermodynamic conditions. This approach benefits from the previous force field developments but unfortunately can generate erroneous results due to excessive ion pairing usually observed in these models when applied at high concentrations.⁵ This behavior is typical for force fields benchmarked against a few standard properties such as the ion–water distances or the Gibbs free energy of hydration, but never rigorously tested and/or optimized for simulations of concentrated solutions. The problem mainly originates from ignoring the electronic polarization that is inherently missing in the most common models, which are nonpolarizable.⁶ The exaggerated electrostatic interactions and unphysical ionic aggregations well below the solubility limits are particularly pronounced for multivalent ions or ions with a high charge density such as lithium or highly charged oxyanions.^7–9 On the other hand, including explicitly polarizable electron clouds is not automatically a panacea, and sometimes polarizable force fields give solution structures of similar quality as nonpolarizable ones.¹⁰ 

In recent years, new ideas and strategies have appeared in the development and application of force fields. One of them is the electronic continuum model (also known as the electronic continuum correction, ECC, or simply charge scaling) proposed as a “lifesaver” for standard nonpolarizable force fields.^11,12 The model suggests the reduction of the ionic charge as a mean-field representation of the electronic polarization of the solvent. This simple and physically justified implementation of electronic screening dramatically improves already existing nonpolarizable models.⁵ For example, the charge scaling prevents excessive salt aggregation,^13,14 allows capturing dynamic properties of kosmotropic salt solutions,¹⁵ and resolves overestimated interactions in biological systems⁵ and complex aqueous solutions.¹⁶ At the same time, charge scaling has raised new questions and challenges, including its applicability to dielectrically inhomogeneous systems^17,18 or systems with external electric fields,^19,20 what is the appropriate degree of charge scaling,^15,20 and whether it is necessary to develop an ECC-compatible water model.^5,20,21 Alternatively, one can benefit from the ECC ideas by implicitly weakening electrostatic interactions while using unscaled (“full,” formal) charges instead of effectively charge-scaled models of solvents and ions.²⁰ The concept of integrating the electronic polarization into already existing models has pushed the quality of MD simulations of ionic solutions to a new, fairly promising level.

Among experimental techniques, neutron diffraction with isotopic substitution (NDIS) has been demonstrated as a powerful method in resolving the complex structure of ionic solutions. This technique is among very few experimental tools directly providing the “true” atomic structure of aqueous salt solutions. The neutron total scattering measurements are based on recording diffraction/diffuse patterns observed when neutrons change their momentum (or wave vector q) after elastically colliding with the atomic nuclei.²² The pairwise radial distribution functions (RDFs) of the constituted elements can be further derived from the collected neutron total scattering signal via a Fourier transformation. With the isotope labeling, measurements can be done on structurally identical solutions but different in the isotope of one of the elements. Since the neutron coherent scattering lengths are different for different isotopes, the direct scattering signal difference between the labeled and unlabeled solutions allows the extraction of the pairwise RDFs only involved for the isotopically substituted atoms.²³ This simplifies deconvoluting the experimental signal, which even in the case of simple electrolyte solutions (those consisting of water oxygens, water hydrogens, cations, and anions, i.e., only four different elements) contains ten pairwise correlations contributing to the total scattering pattern. By labeling one of the elements with the isotope, the isotopic substitution method decreases this number to four components only. Additionally, since the neutrons scatter by interacting with the atomic nucleus, the negative and positive coherent scattering length differences in hydrogen and deuterium (b_H = −3.7046 fm and b_D = 6.671 fm) offer the possibility of contrast matching by isotopes. The so-called null water (ⁿH₂O, an isotopic mixture of heavy and light water with the mole ratio of 0.6407 H₂O:0.3593 D₂O) represents a perfect null scattering solvent media and does not produce contribution from hydrogens since the scatterings from hydrogen and deuterium cancel each other.^16,24 At the same time, the first-order difference analysis, depending on the substituted element, neglects certain “hidden” atom–atom interactions potentially significant in the studied solution. For instance, when an anion (or cation) is labeled by isotopes, signals arising from ion-induced water–water correlation completely cancel out.

NDIS measurements are not only a unique and, to a large extent, self-sufficient experimental technique but can also serve a role of a “policeman” to verify the quality of MD force fields.²⁵ In parallel, optimized MD models can deconvolute the NDIS data providing accurate interpretation and atomistic structural correlations. Together, this twofold approach has led to the development and optimization of several MD models based on both NDIS and ECC principles,^24,26–28 including models for CaCl₂^29,30 and our recent work on KNO₃ and ZnCl₂ aqueous solutions.^16,31 Since the isotopic labeling can be performed on various elements, e.g., metal cations,³² chloride anions,^16,27 nitrogens,^31,33 hydrogens,²⁴ and oxygens,³¹ it would be extremely valuable to verify available NDIS experimental measurements and subsequent MD parameterization by another independent set of NDIS data for the same system. This is especially relevant for concentrated solutions, where ion pairing is usually enhanced, and the solution structure is more complex.

In this work, we investigate 7.3 m CaCl₂ solutions, which appear to be an excellent candidate for such research, particularly due to various reports on concentration-dependent Ca²⁺–Cl⁻ ion pairing.^34–36 We combine neutron diffraction experiments with chloride isotopic substitution (Cl-NDIS) and MD simulations with several recent force field parameterizations. The presented results not only reveal the ion pairing phenomena and thermodynamic properties of CaCl₂ solutions but also disclose advantages (“good”), possible drawbacks (“bad”), and unforeseen dangers (“hidden”) of combining NDIS and MD approaches.

NDIS experimental methodology is similar to the one in our recent work on ZnCl₂.¹⁶ Briefly, we performed Cl-NDIS experiments on a pair of 7.3 m calcium chloride aqueous solutions using null water, i.e., ^natCl/³⁷Cl of CaCl₂–ⁿH₂O. The source of ³⁷Cl was from Na³⁷Cl salt powders produced from the Stable Isotope Production Program at Oak Ridge National Laboratory (ORNL). To extract ³⁷Cl from Na³⁷Cl salts, cation exchange reactions using H⁺ loaded resins (AG 50 W X12, Dowex) were applied to convert dissolved Na³⁷Cl to H³⁷Cl solution. To produce ³⁷Cl-labeled Ca³⁷Cl₂ salt, a few grams of CaCO₃(s) powders (>99.0%, C4830, Sigma-Aldrich) were dissolved in the H³⁷Cl acid solution under ambient conditions. Details about each step described above and additional purification steps via precipitation needed to produce anhydrous Ca³⁷Cl₂ salt powders are given in the supplementary material. Preparation for the target 7.3 m concentrations was followed by dissolving anhydrous Ca³⁷Cl₂ salt powders in null water. Note that the same preparation procedures were applied to ^natCl samples to maintain consistency between ³⁷Cl-labeled and ^natCl-natural samples.

Prepared solutions were sealed into quartz nuclear magnetic resonance (NMR) tubes (507-PP-7QTZ, Wilmad-LabGlass) and loaded in the sample shifter carousel at the Nanoscale Ordered Materials Diffractometer (NOMAD) instrument at the Spallation Neutron Source, ORNL. For each solution, scattering data were collected in 20 frames for a total of 10 h (30 min per frame). The first-order difference signal was obtained by the subtraction between two total structure factors $Fq$ from diffraction data involving isotope substitutions. The generated $ΔFClq$ data were then Fourier-transformed into the real-space to produce the first-order difference RDF $ΔgClr$⁠. More details on data reduction and difference analyses are provided in the supplementary material.

We tested several recent MD force fields for CaCl₂,^{14,29,30,37,38} in conjunction with compatible water models, SPC/E (extended simple point charge)³⁹ or TIP4P/2005.⁴⁰ The ionic parameterizations include models with the scaled charges,^14,29,30 and the latter two were fitted to previously published Ca-NDIS data.³² Throughout the text, MD models are referred to as the first author’s last name and publication year, namely, Deublein2012, Mamatkulov2013, Zeron2019, Kohagen2014, and Martinek2018. We also tested several derivatives of the Martinek2018 model. A detailed summary of all these models can be found in the supplementary material.

For simulations with the SPC/E water model, i.e., for all except Zeron2019, we used a simulation protocol similar to the one in Ref. 30. The cut-off for Coulomb and van der Waals interactions was 12 Å. The long-range electrostatic interactions were treated using the smooth particle mesh Ewald method.⁴¹ The Parrinello-Rahman barostat⁴² with 2 ps coupling time was used to maintain the pressure at 1 bar. The long-range dispersion correction for energy and pressure was applied. The temperature of 298.15 K was maintained by the v-rescale thermostat⁴³ with a coupling time of 2 ps. The geometry of water molecules was constrained using the SETTLE algorithm.⁴⁴ For the Zeron2019 ionic model (also known as the Madrid model) based on the TIP4P/2005 water, the simulation protocol was changed in accord with the original paper, i.e., the Nosé-Hoover thermostat^45,46 with 2 ps coupling time was applied, whereas the cut-off for electrostatic and van der Waals interactions was reduced to 10 Å. All simulations were performed using a Gromacs software, versions 5.1.4 and 2020.x.⁴⁷ 

The simulation ∼40 Å cubic box of 7.3 m CaCl₂ solution contained 1500 water molecules and randomly placed 197 calcium ions and 394 chloride ions. For 4 m CaCl₂ simulations, the numbers of ions were reduced to 108 calcium and 216 chloride ions. Simulations were 100 ns long with a time step of 2 fs. The simulation outputs were calculated from the last 80 ns of these simulations (except the self-diffusion coefficient of water, where the last 20 ns were used for the analysis). The standard error was estimated using the block average over five blocks. More details about MD simulations, including protocols of hydration free energy,⁴⁸ osmotic pressure,⁴⁹ and potential of mean force⁵⁰ calculations, are given in the supplementary material.

The normalized first-order difference signal generated from the Cl-NDIS experiments on a pair (^natCl/³⁷Cl) of 7.3 m CaCl₂ null water solutions can be expressed as

Here, the first-order difference RDF, $ΔgClnormr$⁠, represents a linear combination of the partial RDFs involving Cl atoms, i.e., $gCl−Car$⁠, $gCl−Clr$⁠, $gCl−Owr$⁠, and $gCl−Hwr$⁠. This expression follows previous NDIS works,^16,22,23,31 and the numerical prefactors A_Cl−j are calculated on the basis of solution stoichiometry and tabulated coherent scattering length of each element (see the supplementary material for detailed information). As indicated in Eq. (1), the Cl⁻–O_w correlation signal contributes ∼72% to the total observed intensities, followed by the Cl⁻–Cl⁻ (∼20%), Cl⁻–Ca²⁺ (∼8%), and basically negligible Cl⁻–H_w contribution. The availability of previously published NDIS data with calcium isotopic substitution on 4 m calcium chloride in D₂O water (Ca-NDIS; ^natCa/⁴⁴Ca of CaCl₂–D₂O solutions),³² also used for the MD force field parameterizations,^29,30 allows us to check how the proposed solution structures agree with each other and whether the previously parameterized MD models can adequately reproduce both Ca-NDIS and Cl-NDIS signals.

Figure 1(a) shows the experimentally measured Cl-NDIS $ΔgClnormr$ compared to the signals computed using selected MD ionic parameterizations.^{14,29,30,37,38} The given comparison evidently indicates that the Martinek2018 model excellently captures most of the $ΔgClnormr$ features, particularly the Cl–O_w peak at ∼3.2 Å. Since the overall agreement of MD predictions with Cl-NDIS data primarily depends on the quality of a chloride model, see Eq. (1), it is not surprising that the Martinek2018 model, containing Cl⁻ parameters independently optimized based on Cl-NDIS of LiCl solutions,²⁷ shows best agreement among tested force fields. By adopting the same chloride parameters, we would be able to bring the previous MD results obtained with another chloride model to better agreement with our experimental Cl-NDIS data, and almost to perfect agreement in the dominant Cl—O_w peak, see the supplementary material for an example.

However, another crucial detail of the measured signal is the peak (or its absence) at ∼2.75 Å representing the Ca²⁺–Cl⁻ contact ion pairing. Our experimental Cl-NDIS data suggest very few contact ion pairs present in the solution despite high ionic concentration since there is no obvious shoulder or peak present at this distance, and all we observe at distances < ∼2.85 Å are noises due to Fourier termination artifacts. The quantitative predictions from the experimental data are, therefore, also challenging. Nevertheless, all the models except Zeron2019, explicitly parameterized to reproduce various properties over a wide range of concentrations, clearly overestimate the contact ion pairing compared to Cl-NDIS data. At the same time, Zeron2019 has poor agreement with respect to the rest of the features observed in the experimental Cl-NDIS curve. The absence of significant Ca²⁺–Cl⁻ contact ion pairing signals in the Cl-NDIS data likely indicates that Ca²⁺–Cl⁻ interactions are mainly solvent-shared as can be seen from the corresponding computed RDFs in the supplementary material.

The visual analysis of simulated systems, Fig. 1(b), highlights the differences between Martinek2018 and Zeron2019 models. With Zeron2019, only ∼9% of Ca²⁺ have at least one contact ion pair with Cl⁻, whereas the Martinek2018 model predicts around 70% of Ca²⁺ forming a contact ion pair(s) with chloride anions. The overall tendency in forming contact ion pairs can be also seen from the summary of coordination numbers, Fig. 1(c). Although the Deublein2012 model predicts only 0.5 chlorides per one calcium ion, which is the second smallest coordination number for the compared models, both Cl⁻–O_w and Cl⁻–Ca²⁺ peaks are largely displaced compared to the experimental Cl-NDIS.

The ion–ion and ion–water disbalance of the Martinek2018 model can be explained when analyzing experimental Ca-NDIS data used as a reference during the parameterization of Martinek2018 and Kohagen2014 models. Figure 2 shows computed both Cl-NDIS and Ca-NDIS signals with the Martinek2018 model, with additionally plotted main individual contributions to the total signal. The Ca-NDIS dataset was reproduced in accord with the original studies,^29,30,32 and the total scattering signal can be described as

While the Ca²⁺–Cl⁻ peak is well-defined and separated from the rest of the Cl-NDIS signal, the situation is rather opposite for Ca-NDIS, where the ion pairing signal is “hidden” underneath the Ca²⁺–O_w and Ca²⁺–H_w correlation peaks at ∼2.4 and ∼3.1 Å, respectively. This is due to the use of D₂O water in the previous Ca-NDIS experiments, while null ⁿH₂O water is used in our Cl-NDIS experiments. Additionally, the Ca-NDIS reference measurements were performed with 4 m CaCl₂ solutions, i.e., at a lower concentration with likely even fewer ion pairs present in the solution. These two factors eventually led to the development of the MD models that overpredict ion pairing at higher concentrations.

The next step in the comparison of MD models is how other commonly reported properties beyond the structure are reproducible by selected force fields. The MD data collected in Fig. 3 cover thermodynamic and transport properties of the solution often used as target properties for the force field parameterization, namely, solution density, the Gibbs free energy of hydration, self-diffusion coefficient, and osmotic coefficient.

The overall performance of the selected force fields is manifold. For example, the Zeron2019 model, explicitly parameterized to reproduce the concentration dependence of density, is obviously the best force field in capturing this property at the target 7.3 m concentration. Other models, particularly Martinek2018, are worse and underestimate both experimental and Zeron2019 values. At the same time, all models satisfactory reproduce a commonly used experimental estimate of the Gibbs free energy of hydration.⁵² Martinek2018 and Mamatkulov2013 reproduce this estimate slightly better than others, including the Zeron2019 model, which overestimates the experimental value. Note that the values of the Gibbs free energy of hydration directly extracted from MD (⁠$ΔGMDhydr$⁠) were corrected for the electronic polarization in the case of scaled ionic charges (⁠$ΔGelhydr$⁠), which is necessary for comparison with the experimental reference and full-charge models.^52,55

The most interesting situation is observed for self-diffusion coefficient predictions. Deublein2012 and Mamatkulov2013 models have zero ionic diffusion and significantly underestimate water diffusion at a 7.3 m concentration, which indicates that the solution structure is tightly jammed and almost solidified. This jammed dynamics questions the reliability of the results obtained for these models at high concentrations. The models with scaled charges (Zeron2019, Kohagen2014, and Martinek2018) perform much better, particularly Kohagen2014. Note that the experimental value is given for the 7.1 m concentration, that is, the value should be slightly smaller at 7.3 m.. Zeron2019 underestimates the self-diffusion, which follows from the lower self-diffusion coefficient of the corresponding TIP4P/2005 water model (2.08 10⁻⁹ m²/s vs experimental 2.3 10⁻⁹ m²/s at 298 K).⁵⁶ The self-diffusion coefficient of the SPC/E water model, compatible with other ionic force fields, is instead higher, 2.54 × 10⁻⁹ m²/s.⁵⁶ These values for TIP4P/2005 and SPC/E water models do not account for finite-size effects.⁵⁷ In the supplementary material, however, we show that we detected from no to negligible box-size dependency for highly concentrated 7.3 m CaCl₂ solutions.

Due to the poor performance of the force fields with “full” charges with respect to transport properties, i.e., the solution is almost solidified, the osmotic coefficient is not reported for these models. The concentration dependence of the osmotic coefficient is related to the ion pairing phenomena in the solution, as can be seen by a comparison of Figs. 1(a) and 3(d). The comparison of experimental and calculated osmotic coefficients confirms that Kohagen2014 and Martinek2018 overestimate the ion pairing (the osmotic coefficient is lower than in the experiment), while the Zeron2019 model underestimates it, although the performance is excellent for concentrations up to 5 m. The similarity of the osmotic coefficients for Zeron2019 and the experiment also supports our interpretation of the experimental Cl-NDIS signal, i.e., relatively modest ion pairing.

Overall, the results in Figs. 1 and 3 indicate that Martinek2018 and Zeron2019 models, each with their pros and cons, are the best for predicting physicochemical properties at 7.3 m among those tested. Their somewhat opposite performance emphasizes the fundamental limitations of empirical force fields based on Lennard-Jones potentials and constant partial charges. At this point, it seems impossible to improve either of these models to unite their advantages into a single set of parameters. This is due to the fact that shifting the Cl⁻–O_W peak in Zeron2019 to a larger distance to reach better agreement with Cl-NDIS data leads to exaggerated ion pairing as in Martinek2018, see the supplementary material. Moreover, based on the results presented in the original Zeron2019 parameterization, where models for various monovalent (Li⁺, K⁺, Na⁺, Cl⁻) and divalent (Ca²⁺, Mg²⁺) ions were developed, shorter ion–water distances compared to the experimental reference seem a necessary price to reach a satisfactory performance in reproducing thermodynamic solution properties.¹⁴ However, both Martinek2018 and Zeron2019 models (and basically any) have promising potential, i.e., could be “good,” not necessarily “bad,” when categorized based on their possible applications. On the one hand, the Martinek2018 model, which seems an excellent choice for describing short-range interactions at low ionic concentrations,⁵ would fail in simulations of high ionic strength media due to excessive aggregation. On the other hand, the Zeron2019 model, perfectly suited for modeling complex solutions with various electrolytes under various conditions,⁵⁸ can be unsatisfactory for biological simulations, because the dehydration and tight binding of ions to macromolecules are often a dominant force in, e.g., protein–ion interactions, while the ionic concentrations rarely exceed hundreds of millimoles.

Our study also confirms the necessity of charge scaling or any other implementation of electronic polarizability. Although the importance of incorporating polarizability is becoming widely recognized, we still need to identify proper target properties and relevant conditions when a newly developed force field is aimed to reproduce the experimental data. Some target properties indeed do not necessarily require charge scaling, e.g., the Gibbs free energy of hydration, which can be easily matched using standard nonpolarizable models parameterizing only Lennard-Jones interaction potential and keeping the ionic charge unscaled. This property, together with, e.g., ion–water distances, mostly depends on short-range interactions. That is, correct estimates of these properties can be achieved by adjusting van der Waals forces, regardless of the accuracy of the description of electrostatic interactions. When estimating the Gibbs free energy of hydration, the correction due to electronic polarization must be applied for scaled charges, which unfortunately depends on the choice of an ill-defined ionic radius.^6,59,60 However, the transport properties, such as the self-diffusion coefficient at high concentrations, cannot be adequately resolved without taking into account polarization effects.⁶¹ An increasing ionic concentration also pushes fully nonpolarizable models to their performance limits.

Finally, the value of the charge scaling factor is another debate.^14,15,20 Having an experimental technique as powerful as the NDIS method, it is interesting to assess the correlation between the charge scaling and the ability to reproduce experimental Cl-NDIS measurements usually carried out at high concentrations to increase the intensity of the signal. Here, we explore that possibility using the Martinek2018 model that, at first glance, has the best agreement with both Cl-NDIS and Ca-NDIS data but somewhat suffers from overestimated contact ion pairing. The current version of the model adopts 75% charge scaling in accordance with the principle ECC concept.^6,11,12 The model was optimized to reproduce not only experimental Ca-NDIS data³² but also distance-dependent interaction free energy of Ca²⁺ and Cl⁻ ions in a diluted aqueous solution as benchmarked by ab initio molecular dynamics (AIMD).³⁰ In this work, we prepared several variants of this model applying different charge scaling factors (0.7, 0.8, 0.85, 0.9, 0.95, and 1, i.e., no scaling). The Lennard-Jones potential was also modified in each case so that the positions of Ca²⁺–O_w and Cl⁻–O_w peaks are the same as in the original model.

Figure 4 compares the original Martinek2018 parameterization (denoted as “75%”) with its alternatives in the ability to reproduce both the Cl-NDIS total scattering signal in a 7.3 m solution [Fig. 4(a)] and the AIMD free energy profile of Ca²⁺–Cl⁻ pairing at infinite dilution [Fig. 4(b)]. Apparently, both reference datasets are satisfactorily reproduced by all models, i.e., the models are overall “good.” However, the weakest part remains the overestimated contact ion pairing, despite reasonable agreement of the force field MD free energy with the reference AIMD data. The exaggerated ion pairing can be explained due to the formation of highly coordinated ionic complexes that have more than one chloride in the first solvation shell of a calcium ion, see Fig. 4(c). Recently, we showed that thermodynamic stability and reactivity of ionic complexes depend on the ionic concentration, likely due to changing solvation shells as concentration is increased. This also affects the thermodynamics of every chloride entering the solvation shell of a counter-ion.¹⁶ While AIMD suggests that all Martinek2018-based models adequately describe the binding of the first chloride in a diluted solution, the situation is drastically different for the following chloride binding to the same cation, particularly when the ionic concentration increases. The existence of highly coordinated calcium–chloride ionic complexes in reality is uncertain. Neither Cl-NDIS nor, for example, combined x-ray absorption and x-ray diffraction spectra measurements³⁴ point to significant ion pairing even at concentrations close to the solubility limit. The Zeron2019 model, which seems to capture the ion pairing reasonably well, suggests only ∼8% of CaCl⁺ and ∼1% of CaCl₂ ionic complexes at the 7.3 m concentration. The Martinek2018-like models instead suggest significant populations of CaCl₂, $CaCl3−$⁠, and even $CaCl42−$⁠.

Despite somewhat similar binding of chloride to calcium for all models, in terms of both the ion pairing free energy profile and computed Cl-NDIS signals, the overall hydration of calcium is largely dictated by the charge scaling factor. For example, the number of water molecules in the first hydration shell varies between 6 and 8, when no chloride ion is bound, Fig. 4(d). This difference in the solvation gradually vanishes when one, two, or more chlorides are bound to a calcium ion, and eventually, the coordination number of the first solvation shell converges to a total of 6. Nevertheless, these varying levels of hydration for fully solvated and single-chloride-coordinated calcium evidently affect the structure and behavior of the entire solution, e.g., its density.⁶² 

Given that, we evaluated some of the solution properties for these models, Fig. 5, similarly to our previous comparison. For example, the solution density indeed differs, though poorly described by either model, Fig. 5(a). This again points to the almost impossible reproducibility of both short-range, NDIS-benchmarked interactions and some other thermodynamic properties. Eventually, best agreement is observed for a model with unscaled charges, which, at the same time, is probably the most unrealistic in describing the solution structure under conditions close to the solubility limit. Contrary, the Gibbs free energy of hydration, Fig. 5(b), is reproduced by all models extremely well. Our data also highlight the importance of theoretically calculated, charge-scaling-factor dependent electronic component of the Gibbs free energy of hydration, which is either ignored or simply “hidden” in full-charge nonpolarizable force fields. Eventually, this correction beautifully complements the nuclear part directly calculated from MD simulations. Finally, the self-diffusion coefficient, Fig. 5(c), strongly depends on the charge scaling and is unlikely to be adequately reproduced with unscaled ionic charges. To capture the transport properties of ionic solutions, the charge scaling is highly recommended or even strictly required, although the scaling factor is again a matter of debate.^14,15,20

The neutron total scattering technique is a unique scientific tool that allows one to unravel the structure of ionic solutions at an unprecedently accurate atomic level, especially when coupled to molecular simulations. Here, we carried out neutron diffraction experiments with isotopic substitution of chloride in null water with molecular dynamics simulations on concentrated 7.3 m CaCl₂ solutions. Our Cl-NDIS data estimate the contact ion pairing as insignificant even under conditions close to the solubility limit and suggest rather solvent-shared ion pairing. The comparison of Cl-NDIS and MD data shows that recent MD models parameterized against neutron diffraction with calcium isotopic substitution successfully reproduce a significant fraction of the Cl-NDIS signal and adequately capture the ion hydration structure. However, these models overestimate the Ca²⁺–Cl⁻ contact ion pairing since it was primarily “hidden” in the Ca-NDIS signal due to overlapping with Ca²⁺–O_w and Ca²⁺–H_w correlations in the total scattering. The observed excessive ion pairing is also typical for the models designed for simulations of low ionic strength environments. At the same time, models with good performance at concentrations close to the solubility limit fail to reproduce the NDIS measurements accurately. Solution and ionic properties such as solution density, Gibbs free energy of hydration, self-diffusion coefficients, and osmotic coefficients are reproduced by these models with various successes. The introduction of electronic polarization via scaling down ionic charges, in principle, improves the performance of the models. However, we show that scaling charges should be complemented by selecting appropriate target properties during the force field development, which ultimately defines the accuracy and further success of the model.

Overall, our results clearly indicate that despite enormous progress in the force field development during the past decades, MD models still should be treated cautiously. Even an accurate parameterization against evidently good experimental techniques, such as NDIS, may imply only a certain range of suitable applications for the model. When simulating ionic solutions, one should keep in mind the strong dependence of solution properties on various factors, including concentration. The concentration dependence of the solution structure is challenging to capture due to the inherent limits of nonpolarizable force fields. Although scaling ionic charges helps in resolving some drawbacks and adds another fitting parameter during the force field development, specific limits still exist and can be amplified at high concentrations. This suggests that MD models should not be labeled “good” and “bad” but rather categorized into application categories that are based on their initial design and target properties during the parameterization. The range of relevant concentrations seems to be a good criterium for splitting MD models into certain groups. In this case, NDIS measurements are a perfect tool to verify the performance of MD models since they provide accurate atomic data usually collected for highly concentrated solutions, i.e., under conditions challenging for the majority of the force fields. Moreover, several NDIS datasets may be required depending on the solution composition and substituted element, which would help avoid neglecting the characterization of “hidden” segments of the solution structure.

See the supplementary material for the detailed methodology and additional results of NDIS experiments and MD simulations

D.B. and M.P. were supported by the Ministry of Education, Youth, and Sports of the Czech Republic, Project No. LTAUSA17163. C.T. acknowledges the International Max Planck Research School for Many-Particle Systems in Structured Environments hosted by the Max Planck Institute for the Physics of Complex Systems, Dresden, Germany. H.-W.W., N.R., J.C.N., and A.G.S. were supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division. The research at the Spallation Neutron Source at Oak Ridge National Laboratory was supported by the Scientific User Facilitates Division, BES, DOE. The authors also thank Pavel Jungwirth and Hector Martinez-Seara for fruitful discussions related to this research. The sharing of ab initio molecular dynamics data by Tomáš Martinek is also appreciated.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding authors upon reasonable request.