Density functional theory based molecular dynamics study of solution composition effects on the solvation shell of metal ions

We present an ab initio molecular dynamics study of the alkali metal ions Li + , Na + , K + and Cs + , and of the alkaline earth metal ions Mg 2+ and Ca 2+ in both pure water and electrolyte solutions containing the counterions Cl – and SO 42– . Simulations were conducted using different density functional theory methods (PBE, BLYP and revPBE), with and without the inclusion of dispersion interactions (–D3). Analysis of the ion–water structure and interaction strength, water exchange between the first and second hydration shell, and hydrogen bond network and low-frequency reorientation dynamics around the metal ions have been used to characterise the influence of solution composition on the ionic solvation shell. Counterions affect the properties of the hydration shell not only when they are directly coordinate to the metal ion, but also when they are at the second coordination shell. Chlorine ions reduce the sodium hydration shell and expand the calcium hydration shell by stabilizing under-coordinated hydrated Na(H 2 O) 5+ complexes and over-coordinated Ca(H 2 O) 72+ . The same behaviour is observed in CaSO 4 (aq), where Ca 2+ and SO 42– form almost exclusively solvent-shared ion pairs. Water exchange between the first and second hydration shell around Ca 2+ in CaSO 4 (aq) is drastically decelerated compared with the simulations of the hydrated metal ion (single Ca 2+ , no counterions). Velocity autocorrelation function analysis, used to probe the strength of the local ion-water interaction, shows a smoother decay of Mg 2+ in MgCl 2 (aq), which is a clear indication of a looser inter-hexahedral vibration in the presence of chlorine ions located in the second coordination shell of Mg 2+ . The hydrogen bond statistics and orientational dynamics in the ionic solvation shell show that the influence to the water-water network cannot only be ascribed to the specific cation–water interaction, but also to the subtle interplay between the level of hydration of the ions, and the interactions between ions, especially those of opposite charge. As many reactive processes involving solvated metal ions occur in environments that are far from pure water but rich in ions, this computational study shows how the solution composition can result in significant differences in behaviour and function of the ionic solvation shell.


Introduction
A fundamental understanding of the processes controlling the hydration of ions is of importance to low-temperature geochemistry and biological systems.2][3] In the hydration model proposed by Frank and  Wen to explain deviations of ionic solution behaviour relative to pure water, the insertion of an ion in solution radially aligns water molecules in concentric regions characterised by different levels of interaction between water and ion. 4 The water molecules in the innermost region, also known as the first solvation shell (FSS), are tightly bound to the centric ion and exhibit translational and rotational entropy loss compared with pure water. 5][8][9] Ions in aqueous electrolyte solutions are present as free ions only in dilute solutions.Otherwise, cations and anions form contact or solvent-shared ion pairs. 10The structure and dynamics of the ionic solvation shell can, therefore, be significantly affected by the subtle interplay between the level of hydration of the ions, and the interactions between ions, especially those of opposite charge.In more general terms, the composition of the solution, including the balance between ion-solvent and ion-ion interactions, could result in significant differences in behaviour and function of the FSS.Consequently, an accurate characterization of the effect of solution composition on the structure and dynamics of the ionic solvation shell is important to understand and control reactive processes involving ions in solution.
The MD simulation technique has been extensively used to gain insights into the elusive molecular-level processes controlling the properties of aqueous electrolyte solutions. 11,12However, the force field used to describe the ion-water interaction is crucial to obtain an accurate characterization of structure and dynamics in electrolyte solutions. 13,14An important contribution to this field was the development of ab initio MD (AIMD), where forces are derived from the electronic structure, 15,16 usually in the framework of density functional theory (DFT), 17 which provides the capability of studying non-additivity effects in the dynamics of ions solvation shells.AIMD simulations, using the Car-Parrinello and Born-Oppenheimer schemes, of cations [18][19][20][21][22][23] and anions, 24,25 have been mostly conducted on an isolated ion in pure liquid water, with no counterions in solution.Electrolyte solutions, on the other hand, have been subject to only a few AIMD studies, including the characterization of the dissociation of the NaCl in water, 26 the cooperative ionic effects are of the Na + and Cl À on the hydrogen bonding network, 27 and the influence of a static electric fields. 28ere, we report an AIMD study of the of the alkali metal ions Li + , Na + , K + and Cs + , and of the alkaline earth metal ions Mg 2+ and Ca 2+ , in pure water and electrolyte solutions containing the counterions Cl À and SO 4 2À , in order to characterize composition and concentration effects on the structure and dynamics of the ionic solvation shell.We have chosen the above combination of cations and anions because they are commonly found in nature, including in the composition of the sea and groundwater. 29Salt concentrations ranging from 0.9 to 3.8 mol kg À1 were considered to quantify the effect of solution composition on the strength of ion-water interaction and to characterise the structure and low-frequency dynamics of the FSS.
After having introduced details of the simulations, we next report on the effect of the solution composition on the structure and dynamics of the ionic solvation shell.These results are rationalised in terms of the mechanism of water exchange around cations and the strength of ion-water interaction.
The distribution of hydrogen bonds and the reorientation dynamics of water molecules, in the bulk solution and the first solvation shell, are finally presented.

Computational details
Ab initio (Born-Oppenheimer) MD simulations were conducted with the electronic structure code CP2K/Quickstep code (4.1). 30P2K implements DFT based on a hybrid Gaussian plane wave.The Perdew-Burke-Ernzerhof (PBE) 31 generalized-gradient approximations for the exchange and correlation terms were used, together with the Grimme's -D3 dispersion correction, to provide a more accurate description of the structure of liquid water. 32,33Simulations were also conducted using the Becke-Lee-Yang-Parr (BLYP) functional 34,35 and revised form of the PBE (revPBE), 36 with and without the -D3 term.The Goedecker-Teter-Hutter (GTH) pseudopotentials were used to describe the core-valence interactions. 37All atomic species were represented using a double-zeta valence polarized (DZVP) basis set.The plane wave kinetic energy cut off (E cut ) was set to 1000 Ry.The k-sampling was restricted to the G point of the Brillouin zone.The tolerance for the wave function optimization was set to 10 À6 a.u. that allows for a 1 fs time steps with reasonable energy conservation.AIMD simulations were carried out in the canonical (NVT) ensemble (constant number of particles N, volume V, temperature T) using a Nose ´-Hoover chain thermostat to maintain the average temperature at 300 K. Periodic boundary conditions were applied throughout.

Simulation protocol
First, we have conducted classical MD simulation of 64 water molecules in the isothermal-isobaric (NPT) ensemble (constant number of particles N, pressure P = 1 atm, temperature T = 300 K) for 1 ns to generate an equilibrated aqueous solution.The last configuration was used to generate solutions of hydrated metal (M) ions (single M, no counterions, by replacing one H 2 O with one ion (M: Li + , Na + , K + , Mg 2+ , or Ca 2+ ).Aqueous electrolyte solutions with concentrations ranging from 0.9 to 3.8 mol kg À1 were generated by randomly replacing water molecules with cations, M, and anions, X (Cl À or SO 4

2À
) ions, by making sure that the initial configuration did not contain contact ion pairs.Classical MD simulations used the pairwise Lennard-Jones (LJ) potential model developed by Zeron et al. 38 (denoted as Madrid-2019 forcefield) for Li + , Na + , K + , Mg 2+ , Ca 2+ , Cl À , and SO 4 2À in aqueous solution. 38In Madrid-2019, the water molecules are represented by the TIP4P/2005 model, 39 and the monovalent and divalent ions are modelled using charges of 0.85 and 1.7, respectively (in electron units).This potential model allows a very accurate description of the densities of the solutions up to high concentrations. 38The Madrid-2019 has not been developed for Cs + and to generate equilibrated solutions of CsCl we used the LJ potential model from the Amber database. 40For each system, we conducted to equilibrate the volume of the simulation box followed by AIMD (NVT) simulations has been routinely adopted to study the properties of hydrated ions. 18,20,22,42We did not attempt to run AIMD in the NPT ensemble, which could determine the ''correct'' system density at zero pressure, because depending on the DFT method (functional and dispersion correction) equilibrium density of water can change significantly. 43Instead, we chose to fix the system size to match the one obtained using the classical MD simulation of these systems.NPT trajectories must be long enough to assure statistical convergence, and the convergence of pressure for DFT plane wave MD simulations in the NPT ensemble requires a significantly higher basis set cutoff due to grid sensitivity issues when the volume of the simulation cell can fluctuate. 44,45sults and discussion

Solvation structure
Ion-water radial distribution function.Information regarding the structural properties of the coordination sphere of the alkali metal ions (M = Li + , Na + , K + and Cs + ) and alkaline earth metal ions (M = Mg 2+ and Ca 2+ ) in pure liquid water and in the presence of chloride (Cl À ) and sulphate (SO 4 2À ) counterions, can be deducted from the radial distribution functions (RDFs) of the M-O pairs, g MO (r), in Fig. 1, which represent the probability, relative to a random distribution, of finding a water oxygen atom at a distance r from the metal ion.Comparison of the hydration structure of the cations (isolated ion, no counterions) with the results obtained from other AIMD simulations and experimental measurements (Table S2 of ESI †) shows that the PBE-D3 functional with the hybrid Gaussian (DZVP) plane wave (E cut = 1000 Ry) basis set gives an accurate description of the position and amplitude of the M-O RDF first peak as well as the average number of water molecules in the first hydration shell.For example, our results for Na + and K + are in good agreement with recent AIMD simulations conducted with the strongly constrained and appropriately normed (SCAN) meta-GGA functional for complex condensed phase systems. 19A comparison of the ion-water RDFs obtained from AIMD (PBE-D3) and classical MD (Madrid-2019) are also reported in Fig. S1 of ESI.† Fig. 1 shows that the average metal-water distance of the first coordination shell increases when moving down the alkali and alkali earth group.The lithium, magnesium, calcium and, to a lesser extent, sodium ions are characterized by a more rigid and well defined first solvation structure.As to the hydration structure of K + and Cs + , the intensities of the first peaks are lower and their profiles are much less defined.In aqueous electrolyte solutions, the presence of the counterions in solution increases the intensity and decreases the full-width-half-maximum of the first peak of Li + and Mg 2+ ; no significant effect is observed for the g MO (r) profiles of Na + and Ca 2+ .On the other hand, the chloride ions broaden the RDFs of the K + and Cs + cations, indicating a perturbation of their coordination shells and a more labile hydration structure.
Ion-water coordination shell distribution and ion pairing.We report in Fig. 2 the probability distribution of the number of water molecules in the first coordination shell of the metal ions, obtained from the analysis of the AIMD trajectories of hydrated ions and of the 1.9 kg mol À1 aqueous electrolyte solutions.For the hydrated ions the average size of the FSS, as well as the number of accessible hydration states, increases as we move down the periodic table.The tetra-hydrated lithium complex, [Li(H 2 O) 4 ] + , and the hexahydrated magnesium and calcium complexes, [Mg(H 2 O) 6 ] 2+ and [Ca(H 2 O) 6 ] 2+ , were the only species detected for these ions.On the other hand, hydrated Na + , K + and Cs + can access multiple hydration states; the caesium ion, for example, can coordinate between four and nine water molecules.In aqueous electrolyte solutions, metal ions are not purely hydrated, but they tend to form contact ion pairs (CIPs), solvent-shared ion pairs (SSHIPs), and solvent-separated ion pairs (SSIPs) with the counterions in solution.The percentages of CIPs, SSHIPs and SSIPs in the 1.9 mol kg À1 solutions are listed in Table 1, while the results for the more dilute (0.9 mol kg À1 ) and concentrated (3.8 kg mol À1 ) solutions are reported in Table S3 of ESI.† For the 1.9 mol kg solutions, Cs + has a higher tendency of forming CIPs with Cl À (B73%) than Na + and K + , which prefer solvent-shared or solvated-separated ions pairing.The caesium-chloride RDFs (Fig. S2 of ESI †) shows that the presence of an intense peak centred at 3.6 Å, which is roughly the sum of the ionic radii of Cs + and Cl À in water. 46AIMD simulations of a larger simulation box, containing 8 Cs + and 8 Cl À in 710 water molecules (17 ps), show similar ion pairing (80.0%CIP, 16.7% SSHIP and 3.3% SSIP) and Cs-Cl RDF (3.6 Å) results (Fig. S2 of ESI †).A recent EXAFS study of Cs-Cl ion pairing in dilute and concentrated caesium chloride solutions reported Cs-Cl distances of approximately 3.5 Å, 47 confirming the conclusions obtained in this study.For 1.9 mol kg À1 NaCl(aq) and KCl(aq), the metal-chloride RDFs have a very broad profile between 3 and 4 Å, confirming the tendency to form separated ion pairs.For KCl(aq), CIPs are observed in the more dilute 0.9 mol kg À1 solution (36% CIP).
In the most concentrated solutions, 3.8 mol kg À1 both Na + and K + form CIPs with Cl À , 43.9% in NaCl(aq) and 56% in KCl(aq).The tendency, reported in this study, of Cs + to form CIPs with Cl À contrasts with classical MD simulations of M-Cl (M = Na, K, and Cs) in solution, which reported dissociation rate constants (k) following the trend k NaCl o k KCl o k CsCl . 48One source of this incongruency could be the accuracy of the interatomic potential, since calculations of the dissociation of NaCl in water from AIMD simulations reported a remarkable difference in the free energy profiles between the ab initio and classical potentials. 49,50This stresses the importance of the inclusion of many-body effects in the interaction potentials used to simulate electrolyte solutions. 12Other factors could also affect the population of contact and solvent shared IPs other than the absolute strength of the ion-chloride interaction such as the relative dynamics of exchange of water and chloride around the ions.The analysis of the number of ligand exchange around Na + , K + and Cs + ions shows, in fact, a competitive exchange between the Cl À and H 2 O, whereas for K + and Cs + the frequency of water exchange is significantly faster than chloride (Table S4 in ESI †).
For LiCl(aq), the analysis of the ion-water coordination shell distribution shows that in the presence of Cl À , the number of water molecules coordinated to Li + reduces from four to three as a result of the formation of Li + Á Á ÁCl À contact ion pairs (Table 1).In 1.9 mol kg À1 NaCl(aq), Na + and Cl À form solvent-shared or solvent-separated IP + .Despite not being coordinated to Na + , chloride ions induce a significant contraction of the sodium hydration shell: in NaCl(aq) the sodium ions are in the fivecoordinated state for 90% of the simulation period, compared to approximately 50% of the simulation period of hydrated Na + (isolated ion, no counterions).Consequently, counterions such as Cl À can stabilize under-coordinated metal-water complexes even when they are not part of the first coordination shell of Na + .The Cs + and Cl À ions are, for most of the simulation period, present as CIPs (72%) but due to the high flexibility of caesium first hydration shell, the coordination number distribution is only marginally affected by the presence of chloride ions (Fig. 2).Table 1 The speciation of the cations M (Li + , Na + , K + , Cs + , Mg 2+ and Ca 2+ ) and anions X (Cl À , SO 4

2À
) as a function of concentration: contact ion pair (CIP) when M and X are in direct physical contact; solvent-shared ion pairs (SSHIP) when M and X are separated by one water molecule; solventseparated ion pairs (SSIP) when M and X are separated by at least two water molecules.Assignments made from the analysis of the M-X radial distribution functions (RDF): CIP M-X where r min1 M-X and r min2 M-X are the positions of the first and second minima, respectively, of the M-X RDF.The average CN of the first hydration shell of Ca 2+ is still a controversial issue for AIMD simulations. 51Several DFT based MD simulations, using the Car-Parrinello and Born-Oppenheimer schemes, of the hydrated calcium ion (single Ca 2+ , no counterions) reported six water molecules in the first hydration shell of Ca 2+ , 52 which agrees with our result (Fig. 2).However, experimental results from extended X-ray absorption fine structure (EXAFS) measurements in CaCl 2 solutions, with concentrations ranging from 0.12 to 6 mol kg À1 , resulted in CNs between 6.8 and 8. 53,54 As previously discussed, the presence of counterions (Cl À and SO 4

2À
), which are normally present in experiments, stabilizes the higher coordination states of Ca 2+ , yielding a computed CN of Ca 2+ in CaCl 2 (aq) and CaSO 4 (aq) (6.9) in much better agreement with EXAFS measurements of CaCl 2 solutions.Therefore, one simple reason for the disagreement between theory and experiment could be the composition of the solution used to conduct the simulations.Other factors such as the size of the simulation box, the length of the simulation period, and the choice of the density functional method have also been considered by conducting AIMD simulations of hydrated Ca 2+ (no counterions) in a box containing 63 and 124 water molecules using the PBE, rev-PBE and BLYP functionals with and without the -D3 dispersion correction (Table 2).
The results show that: a simulation period of 50 ps is sufficiently long to get a convergent distribution of the calcium-water coordination shell; Ca 2+ in 63 and 124 water molecules yield CNs of 6.0 and 6.5; BLYP-D3 gives a CN of 6.9 and a coordination distribution that is shifted to higher coordination states compared the average CN obtained with PBE-D3 (6.0) and revPBE-D3 (6.1).Both the system size and the DFT method can, therefore, influence the coordination distribution and average coordination size.This analysis suggests that solution composition, system size, and level of theory should all be considered when comparing AIMD simulations and experimental determinations of ionic coordination numbers.

Dynamics of the ionic solvation shell
6][57][58] We have characterised the dynamics of the first and second hydration shells of the alkali metal ions (Li + , Na + , K + and Cs + ) and alkaline earth metal ions (Mg 2+ and Ca 2+ ) in pure liquid water and electrolyte solutions using the ''direct'' method proposed by Hofer and co-workers. 59his methodology has been previously applied for, among others, the characterization of ligand exchange between coordination shells of hydrated alkaline earth metal ions and their carbonate and bicarbonate complexes, 21,22,60 as well as the quantification of the water exchange frequency around calcium sites in calcite-water interface 61 and hydroxyapatite nanopores. 62n the ''direct'' method, the MD trajectories are analysed for water molecule movements and whenever a water molecule crosses the boundary of the cation coordination shell its path is followed; if its new position outside or inside this shell lasts for more than t* = 0.5 ps then the event is accounted as a real exchange event (N ex ).The choice of t* = 0.5 ps is based on the average lifetime of a hydrogen bond between the solvent molecules, obtained using femtosecond spectroscopy, 63 and is the standard value used in the applications of the direct method; setting t* to 5 ps gives only small differences to the statistics of water exchange (Table S5 of ESI †).From the number of water exchange events, the MRT of the water molecules in the first and second hydration shell of the ion can be computed as MRT = (CN Â t sim )/N ex , where CN is the average number of water molecules in the first or second coordination shell and t sim is the duration time of the simulation. 41A recent comparison by Wolthers and co-workers of the direct and survival function methods to characterize the calcium dehydration frequencies from MD trajectories, 22 concluded that the direct methods could provide shorter mean residence times than SF, but requires shorter simulation times than SF to obtain statistically meaningful output, and is better suited to analyse shorter trajectories generated by ab initio MD simulations.Table 3 reports the values of N ex in the first and second coordination shell of the metal ions, together with the number of exchange events normalised to 10 ps, N ex /10 ps, which represents a measure of the ''lability'' of the hydration shell of the metal.For the alkali metal ions, the first shell solvent exchange dynamics increases substantially on going from Li + (MRT = 11.8 ps) to the other alkali metal ions Na + (1.1 ps), K + (0.4 ps) and Cs + (0.2 ps).The MRT of solvent around H 2 O and Cl À is approximately 0.5 ps, which is comparable to the values obtained for K + and Cs + , confirming the lability of these ions.As reported in Table 3, AIMD simulations of hydrated Mg 2+ (single ion, no counterions) showed no exchanges around Mg 2+ , compared with compared to B50 of such exchanges (more facile dehydration) in 500 ps around aqueous Ca 2+ (61.2 ps).The reactivity of the cations varies, therefore, as follow: Mg 2+ { Ca 2+ o Li + { Na + o K + E Cs + .We have rationalized this behaviour by conducting classical metadynamics MD simulations to compute the free energy as a function of the metal-water coordination number (CN) for Li + , Na + , Mg 2+ and Ca 2+ , from which thermodynamic information on the accessible coordination states for these cations can be extracted (Fig. S3 in ESI †).Ca 2+ , Li + and Na + have several coordination states corresponding to a minimum on the free energy profile and for Na + the transition between these states is almost barrierless.On the other hand, the only coordination state accessible at room temperature to Mg 2+ is six and the free energy barrier to go from the six to the five coordinate state (32 kJ mol À1 ) is about three times higher than the free energy barrier of hydrated Ca 2+ to interchange between the six and seven coordination environments.As Mg 2+ is strongly hydrated, water exchange is drastically retarded in its FSS.
As to the dynamics of the second hydration shell, for all ions, the MRTs is less than 1 ps and close to the MRT of solvent around H 2 O (0.5 ps).This result contrasts with previous values obtained from AIMD simulations of the alkaline earth metal ions Mg 2+ , Ca 2+ and Sr 2+ , which reported MRTs between 3.3 and 6.8 ps.These simulations, however, used smaller simulation boxes than the one used in this study and did not employ dispersion corrected DFT functionals. 21Regarding the effect of counterions in solution to the dynamics of exchange around the cations, in LiCl(aq) the formation of Li + /Cl À CIPs increases significantly the MRT around Li + to 20 ps.The other, very labile, alkali ions Na + , K + and Cs + are substantially unaffected by the presence of counter-ions in solution and the dynamics of water exchange around Na + , K + and Cs + remain very fast.AIMD simulations of MgCl 2 (aq) and MgSO 4 (aq) show no exchanges around Mg 2+ .However, notice the contrasting behaviour of the chloride and sulphate ions on the solvent exchange dynamics around Ca 2+ .In both solutions the counterion form SSHIPs with Ca 2+ , but in CaCl 2 (aq) the MRT is 10.6 ps (faster water dynamics) compared with 61.2 ps (slower water dynamics) in pure liquid water; in CaSO 4 (aq) there were no accounted water exchange events around the FSS of the calcium ions.
A further comment concerns the effect of the level of theory and system size on the solvent exchange dynamics around Na + , K + , Ca 2+ and H 2 O (Table 4).We conducted calculations using different DFT methods, with and without the -D3 correction, and simulation boxes.While the absolute number of exchanges (N ex ) varies significantly, the values of the mean residence time in the solvation shell of Na + (B1 ps), K + (B0.4 ps) and H 2 O in the first and second coordination shells of the metal ions Li + , Na + , K + , Cs + , Mg 2+ and Ca 2+ , with a duration of more than t* = 0.5 ps, obtained from the analysis of AIMD (PBE-D3) simulations of hydrated metal ions (single M, no counterions), aqueous electrolyte solution, pure water and hydrated chlorine ion.Values of N H 2 O ex normalized to the number of cations in solution.Mean residence time (MRT) of water molecules in the first hydration shell computed using the relation MRT = (t sim Â CN)/N ex , where t sim is the period of simulation (in ps) and CN is the average ion-water coordination number of the first and second coordination shell, which was obtained from the integration of the ion-water radial distribution function up to the first and second minimum, respectively.Concentration (b) in mol kg   4).Recent estimates using different interatomic potential models and techniques to characterise the dynamics of water produced give MRT values lying in the range of 60-80 ps, in good agreement with the estimate obtained with the PBE-D3 functional. 22Using the same functional (PBE-D3), the MRT value for Ca 2+ reduces to 13.1 ps from 61.2 ps on increasing the system size, which reveals that for Ca 2+ the MRT can be sensitive to the system size.

Caging effect in the hydration shell of ions
Towards resolving the origin of the different local ion-water interaction, the velocity-autocorrelation functions (VACFs) of the metal ions (M) were computed. 64VACF is defined as follows: where v i is the velocity vector of atom i, N O and N M are the number of time origins spaced by t and number of metal ions respectively.Fig. 3 displays the VACF profiles obtained from AIMD (PBE-D3) simulations of the hydrated cation (single M n+ , no counter-ions) and aqueous electrolyte solutions (1.9 mol kg À1 ).Comparison of the VACF profiles obtained from the full set of NaCl, KCl and CaSO 4 solutions are reported in Fig. S4 of ESI.† MD studies of amorphous systems have used VACF analysis to probe the strength of interaction of network formers (e.g.Si, Al, P) and network modifiers (e.g.Ca, Na) with the surrounding ''cage''. 65,66In particular, the occurrence of a dip to negative values in VACF profiles results from the so-called cage effect for the tagged particle; it takes some time for the particle to escape from the cage formed by its surrounding neighbors. 67The VACF is very sensitive to the local ion hydration environment and changes in the profiles of this time correlation function, such as oscillatory behaviour and position of the first minimum (t 0 ), can be used as descriptors of the strength of ion interaction with the surrounding water molecules and of the rigidity of the interatomic vibration. 68Rigid intra-atomic vibrations lead to a fast-oscillatory trend for the VACFs of Mg 2+ , Li + and to a lesser extent Ca 2+ (Fig. 3).Conversely, Na + , K + and Cs + ions show a much broader first minimum, which is not followed by a marked oscillatory behaviour as their interaction with the surrounding cage is weaker.The VACFs of Na + (t 0 = 0.1 ps), K + (t 0 = 0.16 ps), and Cs + (t 0 = 0.5 ps) in particular, display a smooth decay (Fig. 3A), which suggests a weak caging effect by the solvent; this explains the facile kinetics of dehydration of these cations.The lithium ion shows a much sharper first minimum (t 0 = 0.04 ps) followed by a marked oscillatory behavior because of the more rigid interatomic vibration   between the smaller Li + with the surrounding hydration cage.The comparison of the Ca 2+ and Mg 2+ VACFs (Fig. 3B) shows a very different oscillatory behaviour (t 0 = 0.01 ps for Mg 2+ and t 0 = 0.06 ps for Ca 2+ ), which is indicative of the strengthening of ion-interaction with the surrounding water molecules and increased rigidity of the interatomic vibrations ongoing from Ca 2+ to Mg 2+ .The presence of chloride counterions has no significant effect on the VACF profiles of the alkali metal ions and of the calcium ion.On the other hand, the inset of Fig. 3B shows a smoother decay of the VACF of Mg 2+ in MgCl 2 (aq), which is a clear indication of a looser inter-hexahedral vibration.Therefore, chlorine ions located in the second coordination shell of Mg 2+ weaken the Mg(H 2 O) 6 2+ water 'cage' of the ion.This dynamic flexibility further highlights the effect of counterions on the dynamics of water molecules around Mg 2+ .A final comment concerns the sensitivity of the VACF on the density functional method used in the AIMD simulations.The VACF profiles of hydrated Ca 2+ obtained from AIMD simulations using PBE, revPBE, and BLYP, with and without D3 correction, shows some differences (Fig. S5 in ESI †).The PBE method generates the more marked oscillatory behaviour, to which corresponds a more rigid Ca-water coordination environment compared to the other DFT methods, explaining the absence of water exchange events around Ca 2+ during the AIMD (PBE) simulation (Table 3).

Hydrogen bond statistics
Distribution of hydrogen bonds in bulk solution.The AIMD trajectories of pure liquid water and the aqueous electrolyte solutions were scanned to determine the effect of ions on the hydrogen bond (HB) structure, in terms of the percentage of molecules having n water-water HBs and the average number of HBs (n HB ) per molecule (Table S6 of ESI †).We adopted the following geometrical criteria to determine the existence of a hydrogen bond between two water molecules: 69,70 (i) the oxygen-oxygen distance is less than 3.5 Å; (ii) the intermolecular hydrogen-oxygen distance is less than 2.45 Å; (iii) the oxygen-oxygen-hydrogen angle is less than 30 degrees.The average number of HBs in bulk water computed from AIMD (PBE-D3) simulations is 3.75, which is in the range of experimental values (3.5-3.9)obtained using several techniques. 71In comparison, the revPBE-D3 and BLYP-D3 functionals give n HB values of 3.40 and 3.61, respectively.The average number of HBs reported in Table S6 (ESI †) decreases with the increase of salt concentration, a result in agreement with previous MD studies of aqueous alkali halide salt solutions. 27,60,69However, the extent of the effect on the HB network is specific to the metal ion.In the 1.9 mol kg À1 electrolyte solutions containing the doubly charged Mg 2+ and Ca 2+ ions, the average number of HBs is three, which is significantly lower than the 1.90 mol kg À1 solutions of LiCl (3.53), NaCl (3.40), KCl (3.37) and CsCl (3.50).Fig. 4 compares the percentage of water molecules that engage in n HBs for the solutions considered in this study.In pure liquid water, molecules are mostly engaged in four (70%) and three (20%) HBs, with a similar distribution being observed in LiCl(aq), NaCl(aq), KCl(aq) and CsCl(aq).On the other hand, in MgCl 2 (aq) and CaCl 2 (aq), the distribution of HBs is significantly shifted to lower n, with an almost equal proportion of molecules forming two, three and four HBs.
Hydrogen bond statistics in ionic solvation shell.Gua `rdia et al. revealed significant modifications in the hydrogen bonding architecture for the water the first ionic solvation, which do not persist beyond the second shells. 72To investigate the origin of the ion-specific effects to the hydrogen bond structure, we have computed the distribution of HBs for the water molecules that, at each step of the AIMD simulation, are coordinated to the metal ion.Table 5 compares the percentage of molecules having n water-water HBs ( f n ) and the average number of HBs (n HB ) in the first solvation shell of the monovalent (Li + , K + , Na + , and Cs + ) and the divalent (Mg 2+ and Ca 2+ ) cations, together with the HB statistics for pure water and the water molecules coordinated to Cl À , with the latter obtained from AIMD of hydrated Cl À .For the hydrated metal ions (single M, no counterions), water molecules coordinated Mg 2+ (n HB = 2.20) and Ca 2+ (2.30) display very large deviations from the HB distribution of the bulk (3.75).In comparison, the values of n HB in the FSS of the more labile K + and Cs + are 3.06 and 3.24, respectively.The water molecules coordinated to Cl À (3.61) display only minor deviations from H 2 O. Fig. 5A highlights that the average numbers of HBs in the FSS of metal ions consistently reduces in chloride-bearing solutions compared to hydrated metal ions (single M, no counterions).The interaction between ions of opposite charge generates important differences in the distribution of HBs in the FSS.Speciation analysis of the aqueous electrolyte solutions has shown that metal ions are present as free ions only in very dilute electrolyte solutions. 73Otherwise, ions tend to form contact and solvent shared pairs with the anions present in solution (Table 1).In electrolyte solutions, we can, therefore, consider that water molecules exist in different hydration states, depending on their relative position (coordination) to the cations and anions in solution.Classical MD simulations of MgCl 2 solutions with concentrations ranging from 0.3 mol kg À1 (4 MgCl 2 units in 717 water molecules) to 2.8 mol kg À1 (33 MgCl 2 units in 633 water molecules) were also conducted to determine the fraction of water molecules that are in the bulk (or free water), coordinated to one single ion, or coordinated to both Mg 2+ and Cl À .Fig. 5B shows that about 20% of water molecules are coordinated to both Mg 2+ and Cl À , making cooperative effects important for most concentrations.
Despite the influence to the HB network should be mainly ascribed to the specific Mg 2+ -water interaction, the HB structure of the first hydration shell is governed by the subtle interplay between the level of hydration of the ions, and the interactions between ions, especially those of opposite charge.In more general terms, the cooperative effect of cations and anions in solution reflects a balance between solute-solvent and solute-solute interactions since relatively minor changes to the chemical characteristics result in dramatic changes to the HB distribution.

Water reorientation dynamics
The rotational dynamics of the water dipole, m, was computed from the 1st order Legendre polynomial time correlation: 74 where m(0) and m(t) are the unit vectors defining the orientation of m at times 0 and t, respectively.The averages were computed using several time origins and overlapping intervals, each of length t = 16 ps (Fig. S6-S9 in ESI †).
Fig. 6 reports the concentration-dependent time correlation P 1 (t) profiles of the water molecules, obtained from the AIMD (PBE-D3) simulations of the hydrated cation (single M, no counterions) and of the aqueous electrolyte solutions.The P 1 (t) function starts at 1 and then decays asymptotically to zero because of the random and isotropic reorientation of the water molecules in solution.The early stages of fast loss of correlation are caused by librational motion, whereas the long term decay is due to reorientational motion and can be fit by a bi-exponential model, exp(Àt/t 1 ) + b exp(Àt/t 2 ). 74The overall time associated with this process, t reor , is given by the sum of the fitting parameters t 1 and t 2 .The value of t reor for pure liquid water (64 H 2 O), computed at the PBE-D3 level of theory is 5.3 ps, is in good agreement with the experimental value of 4.8 ps. 75n electrolyte solutions, the presence of solvated ions accelerates the reorientation dynamics of water dipoles (faster P 1 (t) decay) compared to pure liquid water.However, the P 1 (t)  This journal is © the Owner Societies 2020 profiles of NaCl(aq) and KCl(aq) do not consistently decrease with the concentration as they are also influenced by the type of ion pairs present in solution.For NaCl(aq), the fastest reorientation dynamics (fastest P 1 (t) decay) is observed for the 0.9 mol kg À1 solution, which only contains free Na + and Cl À ions.On the other hand, the 1.9 and 3.9 mol kg À1 NaCl solutions, which contains contact and solvent-shared pairs, have a water reorientation behaviour similar to that of hydrated Na + (Fig. 6).The correlation function P 1 (t) is, therefore, sensitive to solution speciation and can provide insights into the combined effect of cations and anions on the water reorientation dynamics.Chowdhuri and Chandra showed that in the vicinity of ions the diffusion coefficients are reduced and the orientational relaxation times are increased compared with pure water. 76We have computed the reorientation dynamics of the water molecules in the first and second coordination shells of the metal ions, and of the water molecules in the bulk (beyond second coordination shells of all ions in solution).Fig. 7 compares the correlation profiles obtained from the analysis of the AIMD trajectories of the hydrated metal ions K + and Mg 2+ and of 1.9 mol kg À1 aqueous KCl(aq) and MgCl 2 (aq).For hydrated K + , the reorientation dynamics of water dipoles in the hydration shells is similar to bulk behaviour.The reorientation dynamics of water dipoles is faster in chloride-containing solutions compared to that observed for the hydrated single metal ion: the P 1 (t) profiles of the water molecules in the first shell of K + in KCl(aq) are drastically accelerated compared to the behaviour around hydrated K + (single ion, no counterions).Conversely, the decay of the P 1 (t) profiles of the water molecules directly coordinated to Mg 2+ are drastically slower than in the bulk, but the effect of solution composition (presence of other Mg 2+ and Cl À in solution) is not significant.This can be explained by the much stronger ion-water interaction of Mg 2+ compared to K + .The correlation function P 1 (t) is, therefore, sensitive to solution speciation and provide insights into the combined effect of cations and anions on the water reorientation dynamics.

Conclusions
We have conducted DFT based molecular dynamics simulations of the alkali metal ions Li + , Na + , K + , Cs + and of the alkaline earth metal ions Mg 2+ and Ca 2+ , in pure water and electrolyte solutions containing Cl À and SO 4 2À , with salt concentrations ranging from 0.9 to 3.8 mol kg À1 , to determine the effect of solution composition on the structural and dynamic properties of the first solvation shell of metal ions.Based on the results we conclude the following: In aqueous electrolyte solutions, the presence of the Cl À and SO 4 2À counterions does not have a significant effect on the ion-water radial distribution functions of profiles of Li + , Na + , Mg 2+ and Ca 2+ .The presence of counterions in solution influences the size of the hydration shell not only because of the formation of contact ion pairs, which reduces the coordination sites available around the metal ion, but also when cations and anions form solvent-shared ion pairs.We found, for example, that in CaCl 2 (aq) and CaSO 4 (aq), the chlorine and sulphate ions in the second coordination shell of Ca 2+ stabilize the seven-coordinated Ca(H 2 O) 7 2+ complex over the six-coordinated Ca(H 2 O) 6 2+ species.
Calculation of the velocity autocorrelation function (VACF) of the metal ions has been used to probe the strength of the local ion-water interaction and the rigidity of the interatomic vibration.We found a smoother decay of the VACF of Mg 2+ in MgCl 2 (aq), which is a clear indication of a looser interhexahedral vibration.Chlorine ions, located in the second  coordination shell of Mg 2+ weaken, the Mg(H 2 O) 6 2+ water 'cage' of the ion.This dynamic flexibility further highlights the effect of counterions on the dynamics of water molecules around ions.Analysis of the hydrogen bond statistics in ionic solvation shell shows how the influence to the water-water network cannot only be ascribed to the specific Mg 2+ -water interaction, but also to the subtle interplay between the level of hydration of the ions, and the interactions between ions, especially those of opposite charge.

Fig. 2
Fig.2Probability distributions of the number of water molecules in the first coordination shells of the metal ions Li + , Na + , K + , Cs + , Mg 2+ and Ca 2+ .Results obtained from the analysis of AIMD (PBE-D3) trajectories of hydrated metal ions (single M, no counterions) and of aqueous electrolyte solution (1.9 mol kg À1 ).

(
B0.3 ps) is only marginally affected by the density functional and number of water molecules in the simulation box.For Ca 2+ , no exchanges were accounted with the PBE functional without the -D3 dispersion correction, and the value of MRT varies from 14 ps (revPBE-D3) to (71.7 AE 16.5) ps (PBE-D3) (Table AE 19.3 b 0.4 AE 0.0 b revPBE-D3 4.1 50 938.7 AE 12.1 b 0.2 AE 0.0 b BLYP-D3 4.5 50 822.4AE 20.6 b 0.3 AE 0.0 b H 2 O 124 PBE-D3 4.2 50 592.7 AE 42.3 b 0.4 AE 0.0 b a The standard deviations of N ex around Ca 2+ have been computed from the variation of overlapping block averages of trajectories lasting 50 ps.b The standard deviations of N ex around H 2 O have been computed from the averages of the number of exchanges around each water molecule in the simulation box.

Fig. 4
Fig. 4 Percentage of water molecules that engage in n hydrogen bonds (HBs) in pure water and electrolyte solutions.Values obtained from AIMD (PBE-D3) trajectories.

Fig. 5 (
Fig. 5 (A) Average number of hydrogen bonds (HBs) in the first hydration shell of metal ions obtained from AIMD (PBE-D3) simulations of hydrated metal ions (single M, no counterions) and of aqueous electrolyte solutions (1.9 mol kg À1 ).(B) Distribution of the number of molecules among the water subpopulations, obtained from classical MD simulations of MgCl 2 solutions.

Fig. 6
Fig. 6 Profiles of the correlation function P 1 (t) for the water molecules obtained from AIMD (PBE-D3) simulations of hydrated metal ions (single M, no counterions) and of aqueous electrolyte solutions at different concentrations.

Fig. 7
Fig.7Profiles of the correlation function P 1 (t) for the water molecules in the first and second shell of K + and Mg 2+ obtained from AIMD (PBE-D3) simulations of hydrated metal ions (single M, no counterions) and of aqueous electrolyte solutions (1.9 mol kg À1 ).
point of 50 ps of AIMD.Details of the solutions considered in this study (number of ions and water molecules, cell lengths of the NVT simulations, solution concentrations) are reported in Table S1 of ESI.† 41 For comparison purposes, we have also conducted AIMD simulations of pure liquid water (64 H 2 O, 1 g cm À3 ) and of Cl À in 63 H 2 O molecules.The procedure of first conducting classical MD (NPT) simulations This journal is © the Owner Societies 2020 Phys.Chem.Chem.Phys., 2020, 22, 16301--16313 | 16303 as the starting Results from AIMD (PBE-D3) simulationsThisIn MgCl 2 (aq) the hexahydrated [Mg(H 2 O) 6 ] 2+ is the only hydrated complex present in solution, and in MgSO 4 (aq) only four water molecules are directly coordinated to Mg 2+ because the magnesium and sulphate ions form very stable contact ion pairs (Table1).Static DFT calculations of the free energy of formation (DG) of the Mg 2+ /SO 4 2À and Mg 2+ /Cl À contact ion journal is © the Owner Societies 2020 Phys.Chem.Chem.Phys., 2020, 22, 16301--16313 | 16305exclusively SSHIPs (99.6%) and the hepta-hydrated calcium complexes are stabilised by the sulphate ions located in the second hydration shell of metal ion.

Table 2
Probability distributions of the number of water molecules in the first coordination shells of hydrated Ca 2+ (single M, no counterions) obtained using different DFT based AIMD methods.N H 2 O /Ca 2+ is the number of water molecules per calcium in the simulation cell; t sim is the simulation period used to compute the coordination number (CN) distribution around Ca 2+

Table 4
Number of accounted water exchange events N H 2 O ex in the first coordination shell of Na + , K + , Ca 2+ , and H 2 O with a duration of more than t* = 0.5 ps.The average coordination number (CN) of first shell is given by the integration of the metal-oxygen radial distribution function up to the first minimum.The values of N H 2 O ex are normalized to the number of metal ions in solution.Mean residence times (MRT) are in ps N H 2 O Method CN t sim N H 2 O

Table 5
The distribution of the number of hydrogen-bonds (HBs) per water molecule in pure water and in the first solvation shell (FSS) of Li + , Na + , K + , Cs + , Mg 2+ and Ca 2+ obtained from AIMD (PBE-D3) simulations of hydrated ions and electrolyte solution.The values are percentages of H 2 O with the given number of HBs.Concentration (b) in mol kg À1