Quantum chemical investigation of the effect of alkali metal ions on the dynamic structure of water in aqueous solutions

We report quantum chemical molecular dynamics (MD) simulations based on the density-functional tight-binding (DFTB) method to investigate the effect of K+, Na+, and Mg2+ ions in aqueous solutions on the static and dynamic structure of bulk water at room temperature and with various concentrations. The DFTB/MD simulations were validated for the description of ion solvation in aqueous ionic solutions by comparing static pair distribution functions (PDFs) as well as the cation solvation shell between experimental and available ab initio DFT data. The effect of the cations on the water structure, as well as relative differences between K+, Na+, and Mg2+ cations, were analyzed in terms of atomically resolved PDFs as well as time-dependent Van Hove correlation functions (VHFs). The investigation of the VHFs reveals that salt ions generally slow down the dynamic decay of the pair correlations in the water solvation sphere, irrespective of the cation size or charge. The analysis of partial metal–oxygen VHFs indicates that there are long-lived correlations between water and Na+ over long distances, in contrast to K+ and Mg2+.


Introduction
Aqueous ionic solutions play an essential role in dynamic chemical processes occurring in nature 1-3 as well as human technologies. [4][5][6] The ions are mobile carriers of electronic charge that can carry out a variety of dynamic functions such as signaling purposes or osmotic pressure changes. The ability to transport electric charges is widely used in chemical energy storage applications such as rechargeable batteries and capacitors, 4-6 with device power and energy density depending greatly on the electrolyte's ion diffusion 7,8 and viscosity 9 properties. To advance the performance of capacitors and batteries, it is therefore crucially important to understand the effect of ionic salt addition on the dynamic structure of water at the microscopic level. However, most of the earlier theoretical [10][11][12][13][14][15] and experimental investigations 11,15-18 of ion-specic properties in aqueous solutions have been limited to self-diffusion properties, static solvation shell structure, e.g. effective coordination numbers (ECNs), 14,15 and pair distribution functions (PDFs). A long-standing discussion in the literature regarding ion-water interactions has focused on the effect of ions on the hydrogen bonding network: 10,11,18 smaller ions with higher charge densities are binding the water molecules in the rst solvation shell more strongly, and thereby enhance the hydrogen bonding network around them (so-called "structure makers" or "cosmotropes"), whereas larger ions with lower charge densities exhibit the opposite behavior (so-called "structure breakers" or "chaotropes"). However, these concepts are only qualitative, and there is no clear denition by which one could quantitatively classify specic ions in these terms. In addition, dynamic properties of the water shell structure around ions has only recently been investigated by time-dependent, correlated waterion spectroscopy and simulations. [19][20][21][22][23] Molecular dynamics (MD) is a powerful computational tool to study the properties of aqueous solutions and provide atomic-scale insights that are helpful for electrolyte analysis. In MD trajectories, the position of every atom is recorded at each point in time, and these time sequences can be used to identify the element-specic physical, chemical, and dynamical properties of the aqueous solutions. 24 However, most MD simulations carried out on aqueous electrolytes are limited in their applicability by the accuracy of the respective classical force elds. 19,23 Certain force elds are incapable of describing the basic features of the liquid and thus fail to reproduce experimentally observed ion-specic correlated motion. 25 A serious drawback of most classical MD force elds is the use of xed atomic point charges to evaluate the Coulomb interaction between atoms, and incorporating more realistic, spatially extended charge distributions with associated polarizabilities involves empirical schemes and incurs signicant computational cost. 26 Polarizable force elds have been developed for the simulation of pure water but have not sufficiently matured to faithfully describe aqueous ionic solutions. As an alternative, so-called "rst principles" MD (FPMD) simulations based on density functional theory (DFT), 27 also oen referred to as ab initio MD (AIMD), 28 can help to overcome the present difficulties of classical force eld methods with ionic charge distributions. Unfortunately, they are computationally rather expensive, 29 making it difficult to obtain trajectories of sufficient simulation time length for relevant model system sizes and complexity. 30 In this work, we instead employ the density-functional tightbinding (DFTB) method, [31][32][33] which is an approximation to traditional DFT that is 100 to 1000 times faster and provides a trade-off between accuracy and computational efficiency. 34 Since DFTB is a quantum chemical method and able to compute nuclear forces during MD simulations by solving Kohn-Sham equations on-the-y, DFTB/MD 35 trajectories not only provide structural features but also contain information on the electronic structure such as charge distributions. This makes the method attractive for the simulation of electrolyte systems. 34,36 The aforementioned pair distribution function (PDF) is an important quantity in the analysis of electrolyte structure. It reports the probability of nding a given particle at a distance R from another particle at R ¼ 0, independent of time. Thus, the PDF can be helpful to understand relative atomic positions and the distribution of the atoms around a particular atom. The PDF can serve as a proper parameter in the solid system where particles do not change their positions with time. However, this "static" PDF does not provide information on the timedependent dynamics of correlations between particles that is reported by the Van Hove correlation function (VHF), 37 a timedependent incarnation of the PDF. It has been shown that the time required for an atom to lose or gain one neighbor is directly linked to the liquid ow. 38 In uids, particles are in constant motion, making local temporal correlations essential to understand the macroscopic transport behavior such as viscosity. For example, the diffusion of ions in aqueous solutions depends on the solution's viscosity, ultimately inuencing the performance of capacitors and batteries that use aqueous electrolytes.
The VHF promises to be a powerful tool in understanding the local and correlated motion of aqueous solution as it relates molecular behavior and the macroscopic transport properties simultaneously as a function of positions (r) and time (t). The decay of the rst peak in VHF displays both ballistic-like motion within the rst solvation shell as well as the correlation between ions. Also, the time decay of the rst peak is linked to the viscosity (h) through Maxwell relaxation time s M . 39 Hence, the VHF essentially links local microscopic behavior i.e., elementary excitation or the time to lose or gain one neighbor (s LC ), to a macroscopic property, namely (s M ). 38,39 Moreover, the VHF can be obtained from inelastic neutron or X-ray scattering experiments via a double Fourier transform of the dynamic structure factor S(Q, u), where Q is the momentum transfer and E ¼ ħu is the energy transfer in scattering. Thus, computation of the VHF provides an opportunity for gaining in-depth molecular-level understanding of element-specic collective motion in an aqueous environment on the picosecond timescale by connecting theoretical simulations with experimental studies. It is non-trivial to resolve element-specic correlation from experiment. 22 To understand the element-specic temporal correlation of water with various cations, in this work, we simulated the VHF in an aqueous solution with different concentrations of chloride salts.
DFTB has already been widely used in the literature to explore properties of pure water, 40,41 biomolecular systems, 42 as well as aqueous electrolytes. 43 To the best of our knowledge however, hydration properties of the aqueous solution of most common salts such as NaCl are yet to be explored with the DFTB approach. Further, the use of non-classical methods in the study of VHF is limited to pure water. 44 Investigations of the VHF of aqueous solution from molecular simulations have thus far been limited to classical MD simulations. 19,25 Here, we performed MD simulations at the DFTB level on the aqueous solutions of the chloride salts KCl, NaCl, and MgCl 2 . Our validation against DFT results shows that DFTB can be used to explore the hydration effect of ions within reasonable agreement with rst principles methods in terms of structure and dynamics predictions. Further, we used DFTB for simulating the Van Hove correlation functions, which show the correlated motion of water molecules and ions. We found that the dynamics of the water hydration shell and the water-cation correlation depends signicantly on the type of salt and their ion concentrations.

Computational methods
DFTB simulations of alkali metal chlorides in aqueous solution were performed using the divide-and-conquer (DC) implementation in the DCDFTBMD soware developed by Nakai et al. 36,[45][46][47] Briey, the DFTB method is based on a two-center approximation of the Kohn-Sham Hamiltonian for valence electrons using optimized electronic parameters as well as interatomic repulsive potentials. The Coulomb interaction between atoms is evaluated using a Taylor expansion of the total DFT electronic energy E in terms of density uctuations dr around the sum of atomic reference densities r 0 . 32 In the DFTB Taylor series expansion, 32 is expressed as: . The series can be truncated at various orders giving different DFTB "avors" with various accuracies in the interatomic coulombic interaction (DFTB1, DFTB2, etc.). In this work, we truncated the DFT energy expression at E 3 , and employed the empirical g-damping for hydrogen bonds, known as third-order DFTB, more conveniently called "DFTB3". 33 The corresponding 3ob parameter set developed by Kubillus et al. 48 was employed in our simulations. The only exception was the repulsive potential for the H-O interaction, where adopted from the 3obw parameter set developed by Goyal et al. 41 In the 3obw parameter set, this particular repulsive potential was specically optimized in an iterative Boltzmann inversion process to reproduce the experimental PDF for bulk water. To further improve noncovalent interactions, as suggested in Goyal's original work, 41 Grimme's D3 empirical dispersion correction with Becke-Johnson damping (D3(BJ)) 49 was employed in our simulations. In the DC-DFTB method, the total system is divided into small, overlapping fragments called subsystems and buffer regions. The subsystem density matrices are computed in parallel using the local DFTB Hamiltonian, and electron population is determined according to the estimated, overall system's Fermi level. Ultimately, the density matrix of the whole system is obtained by summing up the subsystem density matrices. This technique helps to speed up the calculation compared to the standard DFTB method. To facilitate the DC methodology in the DCDFTBMD soware, we selected cubic subsystems with side lengths of 3.0 to 4.0Å and a buffer radius between 5.8 and 6.0Å to ensure the reliability of energetics, as well as the numerical stability of the simulations.
Model systems used in the simulation consist of 1, 2, 3, and 4 ion pairs (cation-anion) with 128 water molecules in a cubic simulation box. The size of the boxes were adjusted to afford the experimental water density under the simulation conditions (0.997 g cm À3 were performed for each model system for a total length of 1.6 ns using the Nosé-Hoover thermostat and a time integration step of 0.5 fs. The rst 100 ps of the trajectory was used for equilibration, and the remaining 1.5 ns of the run was used for analysis, where MD trajectory snapshots were collected every 10 fs. We note that, based on the simulations with different system sizes (1000, 512, 250 and 128 water molecules) for bulk water, Matsumoto et al. 44 found that no signicant system size effects are introduced in the smaller systems, and that 128 water molecules are sufficient to study the VHF and its peak behaviors with time. Also, the typical system size used in the study of aqueous ions with AIMD-based simulations is less than or equal to 128 molecules, with the length of the simulation within a few hundred ps. For example, Galib et al. 15 used 95 water moles for 80 ps, and Wang et al. 14 used 63 water molecules for 50 ps. DC-DFTB scheme used in this study allowed us to achieve 1.6 ns simulations on the system of 128 water molecules with less than 30 K CPU hours.
Here we briey introduce the technique employed to extract the VHF from the MD trajectory. The VHF, G(r, t) between arbitrary particles, atoms or molecules, is dened as: where r is the average number density of particles, N is the number of particles in the system, r i (t) is the position of the ith particle at time t, and d(r) is Dirac's delta function. Hence, the VHF expresses the probability for a particle at the origin r ¼ 0 at time t ¼ 0 nding another particle at distance r, at time t. For the element-specic analysis, one can dene a partial VHF G ab (r, t) between atoms of specic chemical element type a and b as: where N a and N b are the numbers of atoms of species {a} and {b}, respectively. It is worth mentioning that the VHF, G(r, t) in eqn (1) can be decomposed into two quantities: (i) the self-correlation, a measure for particle diffusion, representing the probability of nding a particle i at time t under the assumption that this same particle was at the origin at t ¼ 0, and (ii) the distinct correlation, representing the probability of nding a particle j that is different from i at time t under the assumption that particle i was at the origin at time t ¼ 0, i.e. The latter describes collective motions for distinct particle pairs. The distinct correlation at time t ¼ 0 is identical to the conventional, static PDF, g(r). Van Hove functions were generated and analyzed using the scattering package 50 and functions in the water_-vhf_analysis 51 github repository.

Validation of the DFTB/MD methodology
Since the DFTB method employs parameters obtained by way of tting to selected DFT reference data, it is essential to validate its predictions against known experimental or ab initio results and to gauge the reliability of the chosen DFTB method and parameter sets for the specic materials or chemical systems under consideration. 35 We already noted above that the 3obw parameters for water were shown to accurately reproduce the hydration shell structure and the pair distribution function in agreement with experiment. 41 More recently, the VHF of water was investigated in a study comparing experimental data with different MD simulation methods, including divide-and-conquer (DC) DFTB with the 3obw parameter set and the D3(DJ) dispersion. In that study, the DC-DFTB3-D3(BJ)/3ob method was found to reproduce the overall nature of the VHF and decay characteristics of various peaks in very good agreement with experiment. 44 For the purpose of this work, we rst compared the hydration shell structure of cations obtained from DFTB/MD simulations with the ab initio MD (AIMD) results reported by Wang et al. 14 as shown in Fig. 1. We found generally very good agreement between the DFTB/MD and AIMD PDFs for the position of the cation's rst solvation shells, with differences not exceeding 0.2 A. In the case of the Mg-O PDF, the peak positions are in nearly perfect agreement, and only the intensity of the second peak is somewhat underestimated in DFTB/MD. For both Na + and K + , the height of the rst peak in the DFTB/MD cation-water PDF is likewise somewhat underestimated, whereas here the intensity of the broader second peak is in excellent agreement between DFTB/MD and AIMD. We believe this level of agreement is acceptable since slight deviations in peak position and the height of the PDF are not expected to limit the understanding of the overall trends in the hydration structure of different cations, their effects on the water structure, or the effect of concentration changes. It is further important to recognize that AIMD performed with various DFT functionals can also lead to different PDF peak positions and intensities. 14,52 3.2 Effect of ions on the static structure of water and solvation structure of the ions In our study, we have used chlorine as a counter anion for all three cations in an aqueous solution. At the low salt concentrations employed here, cations and anions are expected to be solvated by water molecules instead of forming cation-anion pairs. We have plotted the Cl-water PDF g Cl-O and its running coordination number (n Cl-O ) in Fig. 2. All cations have very similar n Cl-O with the rst peak located z3.1Å. Only a slight difference is observed when comparing g Cl-O of the MgCl 2 system at a distance between 5.5-6Å to those of the KCl and NaCl systems. This observed feature in MgCl 2 agrees qualitatively with the study performed with Xrays and classical MD simulations of MgCl 2 solutions. 16 Also, the running coordination number of Cl with O (n Cl-O ) in the inset in Fig. 2 is identical for all cations, and it lacks any plateau indicating the continuous distribution of water molecules around chloride ion without dened solvation shell. The similar nature of chloride-water PDF and coordination distributions in three cations suggests that the effect of Cl À in water structure will be the same, irrespective of the cation. We therefore conclude that the specic trend observed among the cations regarding water structure and dynamics pertains to cationic features.
To understand the effect of various cations as well as their concentrations on the static structure of pure water, we have plotted its crucial O-O PDF, g O-O , as well as the O-H PDF, g O-H , in Fig. 3 for all simulated systems. In comparison to g O-H , which seems only affected by salt concentrations in the height of its rst peak at 1.9Å, the g O-O of the ionic solutions shows a more signicant deviation from pure water. Although the position of the rst peak corresponding to nearest-neighbor O-O correlations at 2.8Å remains unchanged with the addition of different salts, the intensity of this peak signicantly decreases with an increase in salt concentrations. This decrease is consistent with experimental observations and was explained previously as a result of the conuence of multiple ion-O pair positions whose positions coincide with the water O-O peak. 23 The greatest overall variation of the O-O PDF as function of salt concentration is observed for K + , with Na + a close second, and the least variation is visible in that for Mg 2+ , which is consistent with the order of the respective ionic radii: the ionic radius of K + is 1.38Å and that of Na + is 0.95Å while it is smallest for the Mg 2+ dication with 0.65Å. Such a trend was previously noted and is consistent with greater solvent exchange reactivity for larger cations. 14 In regards to salt concentration changes, we note that deviations of the g O-O from pure water increases with higher number of ion pairs in the solution, and that these deviations tend to "smear out" the pronounced static water structure with increasing salt concentrations.
The observations above suggest that the salt ion pairs are partially breaking the hydrogen bond structure of water in their vicinity, with Mg 2+ having a less disruptive effect than Na + and K + . In the following, we analyze in greater detail the solvation structure of different metal (M) cations using the cation-water PDFs, g M-O , and the cation-oxygen coordination number n M . The PDFs shown in Fig. 4 represent the distribution of water oxygen atoms at various distances around the cations. The positions of the rst peak in the PDF provides an average of the shortest cation-oxygen distance, r max , which are consistent with the ordering of the cationic radii mentioned above. Among the three cations, the rst hydration shell of K + is furthest away from the ionic nucleus and considerably less dened, with the rst peak in the PDF being broader compared to Na + and Mg 2+ . Correspondingly, the heights of the rst peaks in the g M-O PDFs increases in the order K + < Na + < Mg 2+ . We do not observe a signicant variation of these PDFs with salt concentrations within the range investigated here. Our ndings are in very good agreement with the earlier AIMD results reported by Wang et al. 14 The area underneath the cation-water PDFs, g M-O , when integrated from the cation nucleus outwards up to a given point, represents the running coordination number, as shown in the insets of Fig. 4. In these plots we can clearly observe that, different from Na + and Mg 2+ , the running coordination number of K + increases continuously without exhibiting a well-dened plateau. Such plateaus are indicative of the existence of a stable solvent coordination shell, here the so-called cation hydration shell, and the integrated value of n M-O at the plateau  indicates the number of water molecules that comprise the rst hydration shell. The Mg-O running coordination number reaches a well-dened plateau at 2.5Å indicating the rst hydration shell, and the number of water molecules in its rst hydration shell is 6. The nature of the Na + running coordination curve lies between that of K + and Mg 2+ , but the number of water molecules in its rst hydration shell is also 6. The size of the hydration shell around cations is typically a function of their ionic radii, where smaller ions tend to feature a more well-dened hydration shell and larger ions are associated with more poorly dened hydration shells.
For larger ions such as K + , when there is no clear plateau in the running coordination number, the cation-O coordination number is not sharply dened, which does not allow a clear denition of the number of water molecules participating in the rst hydration shell. For such poorly dened coordination shells, various empirical techniques have been devised to estimate an approximate cation coordination number. 14,15,53,54 In our study we have employed the technique proposed by Hoppe 55,56 who suggested to estimate an average M-O distance and derive from it the so-called effective coordination number (ECN). The average distance, l av , is dened as: where the summation runs over all oxygen neighbors of a particular cation, and l min is the smallest cation-O distance in the i th cation-oxygen polyhedron. The ECN is dened as: Based on this denition of l av , the ECN derived from it represents the ion coordination, based on the weighted average distance of the central cation to its nearest neighboring O atoms. Table 1 shows that the average cation-O distance is largest for K + , decreases for Na + and becomes smallest in the case of Mg 2+ , consistent with their cationic radii. In each case the distances, r max , obtained from the PDFs from the centroid of a tted Gaussian distribution function, are slightly larger than the ones estimated from eqn (3) dening the l av . The deviation between r max and l av is highest for K + , which could be expected based on its more diffuse hydration shell. There is no signicant variation between the cation-O distances with changes in salt concentration, as previously noted for the M-O PDFs. The calculated average distances l av are 2.90Å for K-O, 2.35Å for Na-O, and 2.09Å for Mg-O. We nd good agreement between average distances obtained from our DFTB/MD simulations with previously reported AIMD and experimental results for Na-O and Mg-O, see Table 1. The K-O distance is slightly (z0.1Å) overestimated compared to the experimental value, nevertheless it lies within the range reported from a previous AIMD simulation. Among the three cations considered, the ECN is highest for K + , followed by Mg 2+ which is close to Na + , Table 2. Although the ECNs do not change much with the salt concentration, we nd that alkali metal ECNs are slightly higher at higher concentration (1.6 mol kg À1 ) compared to the lower concentration at 0.80 mol kg À1 . Only in the case of Mg 2+ , its ECN remains almost unchanged z5.64Å within the considered concentration range. Our calculated ECNs at 1.6 mol kg À1 concentration 6.86 for K + , 5.54 for Na + , and 5.64 for Mg 2+ , are in reasonable agreement with values reported in the literature from AIMD and experiments, see Table 2. In the case of Mg 2+ , its ECN is slightly underestimated by z0.36 compared to AIMD and experimental values. On the other hand, its running coordination number at the plateau in Fig. 1(c) indicates a hydration number of 6, which exactly matches the literature values, Table  2. This indicates that Hoppe's denition of an ECN may not give consistent results across ions with different cationic charges.
In summary, we nd that DFTB/MD simulations paint a picture consistent with experimental and available AIMD data, where the largest cation under investigation here, K + , features Table 1 Cation-oxygen distances for the first solvation shell, inÅ, derived from the RDF from the centroid of a fitted Gaussian distribution function, r max , and average distance, l av , calculated according to eqn (3). In addition, available DFT and experimental data (Expt.) are listed for comparison. In their experimental study Mahler et al. 17 and Callahan et al. 16 have used X-ray diffraction and Glezakou et al. 57 and Galib et al. 15  the most exible solvation shell and the smaller Na + and Mg 2+ ions have more rigid solvation shells. On the other hand, the loss of a pronounced water shell structure in the static PDF with increasing salt concentration, is more visible in the case of the alkali metal ions, and least noticeable for the Mg 2+ cation. Nevertheless, our analysis of the effect of salt ions on the static water structure is consistent with that of a structure breaker or chaotrope, in the approximate order K + > Na + > Mg 2+ .
3.3 Effect of ions on the dynamic water structure - VHFs   Fig. 5 shows the total VHF, G(r, t) of pure water and the aqueous ionic solutions of KCl, NaCl, and MgCl 2 salts at a concentration of 1.2 mol kg À1 molality. Notably, a comparison of the water VHF at room temperature between experiment and various simulations on different potentials previously showed that DFTB/MD could reproduce experiment very well. The rst peak of the pure water VHF at t ¼ 0 is at z2.8Å and the second peak is at z4.5Å, consistent with the static PDF of pure water. As time elapses, the rst peak position moves towards higher values of r and the second peak position moves towards lower values of r. 44 As discussed above for the static PDFs, the addition of ions can signcantly inuence the water VHF and even introduce distinct new spectral features. Unless otherwise noted, we focus on a salt concentration of 1.2 mol kg À1 molality for clarity, and is sufficient to illustrate the important effects of the salt ions on the water VHF. As was the case for the static water-water PDFs for different salt concentrations shown in Fig. 3, the total VHFs depicted in Fig. 5  Other changes in the VHF are somewhat more subtle and are better discussed based on the heatmaps in Fig. 5(b). At rst glance it is immediately apparent that the pure water VHF is the spectrum with the most pronounced features, in particular for longer periods of time. This is particularly true for the intensity and long-time decay of the features corresponding to the second and third peaks at z4.3 and 6.5Å. These peaks can be attributed mainly to the O-O pair correlation with molecules in the second and third self-solvation shell. They are less pronounced for the KCl and NaCl salts, with MgCl 2 being somewhat in between the other salts and pure water in terms of these longtime decay features. Second in prominence is a feature that is easily overlooked in the time series plots but clearly visible in the heat maps, namely the delayed decay of the rst peak corresponding to nearest-neighbor O-O correlations at 2.8Å for all salt solution VHF spectra. At time scales above 0.5 ps, this peak merges with the peak intensity originating from intramolecular self-and distinct pair correlations, and this merger is a characteristic feature for all ionic solutions. We further note in regards to the rst peak that over time it seems to develop a slight shoulder structure at z3.21Å which corresponds to interactions between water oxygen and the chloride anion. We note that qualitatively similar observations were previously made in experiment and classical MD simulations for Na + salts with different anions. 22 An interesting observation can be made for Mg 2+ , where this delayed decay of the rst peak is highly prominent, yet the decrease of intensity and long-time decay features for the second and third peaks at z4.3 and 6.5Å are less pronounced.
More information can be extracted from the decay behavior of the rst-neighbor VHF peak, 19 which is related to transport properties and viscosity as they are directly linked to the correlated motions between neighboring particles. 20,21,25 The height, G(t) and area under the curve A(t) of the rst peak is shown in Fig. 6 for pure water and the three salt solutions at  1 are the local minima on the le and right to R 1 , respectively. The decaying behavior can be tted by using a two-step relaxation function. To this end, following the approach described in ref. 22,44, we have tted A(t) (t # 0.75 ps) to a compressed exponential function of the form A(t) ¼ a 1 e (Àt/s 1 )g 1 +a 2 e (Àt/s 2 )g 2 , where s 1 and s 2 represent the relaxation times of the rst-step and second stepdecay. The tting results are listed in Table 3. Unsurprisingly, the values obtained here for pure water are in good agreement with the DFTB/3obw values reported previously. 44 It is notable, however, that our results for the aqueous NaCl solution also agrees with the corresponding experimental results previously reported by Shinohara et al. 23 The addition of NaCl leads to a decrease in a 1 and an increase in a 2 . Likewise, a 1 decreases with the addition of KCl and MgCl 2 , while a 2 decreases for KCl and remains unchanged for MgCl 2 . s 1 remains almost unchanged with the addition of salts, with some decrease associated with s 1 for NaCl. This nding indicates that the rst step decay of pair correlation is not signicantly changed with the addition of salts. However, s 2 increases signicantly with the addition of salts, indicating that the addition of salts induces a slower second step-decay behavior. KCl and MgCl 2 exhibit a similar decay behavior among the three salts, in and order of decay time s 2 decays as Na + z K + < Mg 2+ . An increase in salt concentration will further increase the decay time of the rst peak height, as shown in Fig. 7.
Since the concentrations of the aqueous ionic solutions under study are relatively low, most of the features of the total VHF are dominated by O-O correlations. To analyze the ion-specic correlation function in greater detail, we have also computed partial VHFs and their intensity maps for water-ions, shown in Fig. 8. O-Cl À partial VHFs, Fig. 8(a) and (b) show that the rst peak located between 3 and 4Å is similar in all three cations, which is also the case in the static PDF shown in Fig. 2. Additionally, it is worth pointing out that the O-Cl correlation in the rst peak Fig. 8(b) is extended in time compared to O-O peak correlations observed in pure water VHF Fig. 5(b). Our nding of a slow decay of the water-Cl À rst peak correlation agrees qualitatively with the experimental study performed on NaCl solutions. 23 The second peak between 4 and 6Å shows some ion-specic features. Namely, (i) in MgCl 2 , the second peak is more diffuse in space compared to NaCl and KCl. (ii) Time correlation of the second peak is shortest for KCl, followed by NaCl and MgCl 2 .
Next, we discuss the feature of water-cation VHFs depicted in Fig. 8(c) and (d). All three cations give rise to their rst peak at the sum of ionic radii with the oxygen, namely, 2.78Å for K, 2.35 A for Na, and 2.05Å for Mg. The decay of the rst O-M peak is remarkably slow, as it remains clearly visible at time t ¼ 1.5 ps irrespective of the salt type, indicating a strong association between the cation and the water molecules in the rst solvation shell. The second peak, located between 4 and 6Å, shied towards higher values of r and is short-lived in K + compared to Na + and Mg 2 . The third peak remains only visible for NaCl, which is located at a slightly shorter distance than it appeared in the total VHF of pure water, Fig. 5(d). The partial O-Na VHF indicates that there are long-lived correlations between water and Na + over long distances, unlike in the case of K + and Mg 2+ .
To summarize our observations on the effect of ions on the dynamic structure of water, we nd a common picture Fig. 6 (a) Height of first peak, G 1 (t) as a function of time, (b) area under the first peak A(t) of pure water compared with aqueous solutions of KCl, NaCl, and MgCl 2 at a concentration of 1.2 mol kg À1 molality.

Conclusion
In summary, we have performed DFTB-based MD simulations for aqueous solutions of KCl, NaCl, and MgCl 2 at different salt concentrations. Comparison with available AIMD simulations and experimental results conrmed that DFTB/MD simulations are in the position to accurately describe the effect of cations on the static water structure, and reproduces the AIMD hydration shell of cations in aqueous electrolytes. Our simulations demonstrate for the rst time that DFTB/MD methodologies are allowing accurate simulations at a rst principles level without the need to empirically t charges. This opens the possibility to employ the method to fundamentally understand the local correlated motion of water with salt ions through the VHF dynamic pair correlation function, which is not easily accessible to AIMD simulations due to the required simulation system sizes and timescales. Considering only the static PDFs, it is found that the largest cation under investigation here, K + , features the most exible solvation shell and the smaller Na + and Mg 2+ ions have more rigid solvation shells. The Mg 2+ appears to have the smallest effect on the static PDF of pure water, and we nd the following order of the cations with respect to their ability to disrupt the water network: K + > Na + > Mg 2+ . This is consistent with the fact that pair correlations in the VHFs of the salt solutions are decaying slower than for pure water, irrespective of the type of salt, with Mg 2+ showing the slowest correlation decay. Partial VHF analysis indicates that the O-M correlations within the rst solvation shell of the cations decay remarkably slowly, with the decay of the second solvation shell being least pronounced for Mg 2+ and Na + . Surprisingly, the sodium cation exhibits the longest-lived correlations with water motions around it, indicating more pronounced long-range effects on water dynamics.
Our combined quantum chemical simulation of the effect of cations on the static and dynamic structure of water further gives credibility to the thesis that it is difficult to draw a clear line between chaotrope and cosmotrope ions, as all ions are delaying correlation decay in their solutions, with Mg 2+ having the least impact on the static water structure but greatest impact on its dynamic correlations.

Conflicts of interest
The authors have no conicts to disclose.