Communication: Modeling of concentration dependent water diffusivity in ionic solutions: Role of intermolecular charge transfer

The properties of water and aqueous ionic solutions are of great scientific interest because water plays a central role in the atmosphere, biological environments, and various industrial processes.^1–9 Of particular interest are solutions of relatively simple salts consisting of alkali cations and halogen anions. The translational diffusivity of water in solutions of alkali halide salts depends on the identity of ions, exhibiting dramatically different behavior even in solutions of similar salts of NaCl and KCl.^10,11 While water translational diffusion coefficient in NaCl decreases with the increase of salt concentration, it slightly increases as the concentration of KCl increases. Given the similarity of these two salts, a simple model is unlikely to explain this qualitative difference in the behavior, thus producing a challenge to molecular dynamics (MD) simulations that use classical force fields.¹¹ 

MD simulations that use force fields with permanent charge values (PQFF) have been quite successful in studying homogeneous water and also provided reasonable values for thermodynamic quantities such as energies of ion solvation.¹² For description of water in inhomogeneous environments, such as interfaces and clusters, it is important to include fluctuating character of charges due to electron polarization,^13–16 and this was done using force fields that employed flexible charges (FQFF).^17–20 However, simulations that use either PQFF or FQFF are not able to reproduce the experimentally observed concentration dependent water diffusivity in solutions of alkali-halide salts.¹¹ Recently, Ding et al. showed that First Principles MD (FPMD), which takes dynamical changes of electronic structure into account, is able to describe the different trend displayed in water diffusivity of 3M solutions of NaCl and CsI, while simulations using traditional PQFF and FQFF are not able to do so.²¹ Nevertheless, it is still not clear if classical force field can be modified to reproduce the experimentally observed concentration dependent water diffusivity for solutions of alkali halides and retain discrete representation of electron charges.

Our recent simulations of systems containing water and just one ion have shown that an explicit treatment of charge transfer via FQ-DCT model^22–25 in classical MD simulations yields results for water diffusivity that are in agreement with FPMD simulations.²⁶ We use this model here to investigate the experimentally observed concentration dependent water diffusivity in ionic solutions of NaCl and KCl, and we demonstrate that the MD simulation that uses force field that explicitly considers charge transfer can reproduce the experiments, as shown in Figure 1.

We will start by discussing the extent to which classical MD simulations that use the FQ-DCT (Fluctuating Charge–Discrete Charge Transfer) model can reproduce the essential physical features by comparing to FPMD simulations using the same exchange-correlation (XC) potential (revPBE-D3) and other computational parameters as in Ref. 21. Let us first consider the issue of ion pairing, the propensity of ions to pair: it can be understood by analyzing the behavior of the potential of mean force (PMF) as a function of the cation-anion separation distance. Figure 2 shows that PMF curves obtained from MD simulations using classical PQFF (SPCE-HMN^12,27) model and from FPMD simulations look quite different for both NaCl and KCl. The PMF obtained using PQFF displays deep free energy minima for contact NaCl and KCl pairs, and there exist large free energy barriers separating contact pairs and solvent separated pairs.²⁸ The deep minima and high barriers are not present in the PMFs from FPMD simulations. At the same time, PMF curves calculated using FQ-DCT model reproduce many of the qualitative features obtained from the FPMD simulations. Especially prominent is the reduction of the minimum in the PMF at the ion contact. This may be explained by the reduction of the charge on the ions due to charge transfer and therefore the reduction of the direct Coulomb interaction between the ions when FQ-DCT model is used. Indeed, as Figures 3(a) and 3(b) show, the absolute values of charges on the ions are reduced in simulations using FQ-DCT model and in FPMD simulations. Notice that in both NaCl and KCl cases, the cation charges remain nearly constant as a function of the cation-anion distance. Interestingly, the amount of charge on the anion is largely independent of the identity of the cation in both FQ-DCT model and FPMD indicating that charge transfer between cation and anion does not dictate the amount of charge on the ions. This behavior is rather insensitive to the choice of exchange-correlation approximation we used in our FPMD simulations as shown in Table S2 of the supplementary material.⁵¹ 

Correctly describing charge transfer from the ions to water molecules is essential to simulating the concentration-dependent diffusivity change. As discussed in our previous publication,²⁶ the charge on water molecules plays a significant role on the water diffusivity. Increased charge on water reduces the water diffusivity while decreased charge enhances the water diffusivity. Figures 3(c) and 3(d) show 2-dimensional maps of charge transfer from ions to water molecules. The maps are shown for three different cation-anion separation distances corresponding approximately to the first minimum (contact ion pair), first maximum, and second minimum (solvent separated ion pair) in the PMFs obtained from the FPMD simulations. Figures 3(c) and 3(d) show the charge distribution around NaCl and KCl, respectively. Around the anion, the first shell of water is expectedly charged negative in both NaCl and KCl solutions, indicating that negative charge is transferred from Cl⁻ to water. When the ion pairs are in contact or are solvent separated, the water molecules around Na⁺ ion and K⁺ ion are also negatively charged in the first shell. This is contrary to intuitive belief that main charge transfer occurs between the ion and first shell water, which would produce positive charges on the water molecules. This observation, from both FPMD simulations and the MD simulations with FQ-DCT model, shows that charge transfer is also present between water molecules in the first shell and those farther out, and it clearly evidences that the influence of the ion is not restricted to just the first shell of water molecules. The same observation about charge-transfer from the 2nd hydration shell of the cation overcompensating the charge transfer from the first shell was made in Ref. 23. This shows that presence of the ion is reflected on the charge distribution among water molecules beyond the ion’s immediate vicinity through non-local quantum mechanical effect of electrons.

As we can see, the FQ-DCT model is able to reproduce essential physical features in classical MD simulations as well as FPMD simulations do. Moreover, from Figure 1, we observe that the results from simulations using FQ-DCT model are in a very good agreement with experiment for both NaCl and KCl, while PQFF and FQFF models do not perform satisfactorily.¹¹ Our previous work on simulations containing just one ion solvated by water showed that diffusion of the first shell of water around Na⁺ slowed down, accelerated around Cl⁻, and nearly did not change around K⁺ ion, when compared to bulk water.²⁶ Figure 4 shows the behavior of water translational diffusion in different spatial regions of solutions at different ion concentrations. The water in the first shell of Na⁺ cation slows down substantially for every salt concentration. The figure also shows that water in the first shell of Cl⁻ anion in the 1M NaCl solution is diffusing as fast as in the bulk, but it slows down when salt concentration increases. While in the NaCl salt solution large changes in water mobility occur in all regions defined with respect to the ions, the situation is quite different when the salt is KCl. As seen in Figure 4, for the case of KCl solution, only small changes in water mobility occur in all regions when salt concentration changes. We also observe that diffusion of water around K⁺ ion is slightly slower than in bulk, and it remains rather independent of the salt concentration. These results show that such a qualitatively different behavior in the concentration dependence of water diffusivity in solutions of NaCl and KCl is not a direct result of a difference between the cations alone, since the water diffusivity also changes around Cl⁻ with the salt concentration change, as evidenced here.

It was recently argued that an efficient way to take into account the effect of electronic polarization is to screen the ion-ion interaction by reducing the values of the ionic charges.^29–33 This electronic continuum correction with rescaling (ECCR) model was also used to explain the behavior of water diffusion in aqueous ionic solutions as a function of the salt concentration. When it is used with some specific water model, one can qualitatively predict the difference in the behavior of solutions containing salts such as NaCl or CsI.²⁹ The ECCR model is an effective mean-field model that treats charges on different ions in the same way. At the same time, from our simulations and FPMD, we observe that significant charges are transferred to water molecules rather extensively, and that ionic charges are not of equal value. The charge transfer model is able to take this charge redistribution into account and therefore can tackle such delicate issues as the difference between Na⁺ and K⁺ ions in their influence on water diffusion. The distinct difference for the concentration dependent water diffusivity between NaCl and KCl solutions reflects the importance of many-body effects among the water molecules that are outside the direct contact with the ions and not only those waters that are in the first shell around the ions.

Ding et al.²¹ suggested that treatment of simple electrolyte systems like alkali halides requires explicit inclusion of electronic degrees of freedom within FPMD approach. Our present simulations demonstrate that inclusion of intramolecular charge transfer into classical force field can quantitatively reproduce the difference in the behavior of water diffusion in electrolytes such as NaCl and KCl, thus providing a partial answer to the challenge issued in Ref. 11. The explicit inclusion of charge transfer allows us to capture the highly non-local effect even in classical MD simulations for which the electron charges are discrete. Still, more work is needed to examine the possibility to reproduce experimentally observed water diffusion in alkali halide solutions of different salts by using charge transfer models like FQ-DCT. Also, it will be interesting to study how the implementation of FQ-DCT model correlates with data on the rotational diffusion of water around ions.³⁴ Nevertheless, our results inspire an optimism that charge transfer models may correctly reproduce dynamic heterogeneity of water responsible for the subtle behavior of water diffusion in salt solutions.¹¹ Moreover, availability of classical force fields that correctly reproduce dynamic heterogeneity of water in systems containing large molecules with charges, such as proteins and/or DNA fragments, will make simulations of these systems much more computationally feasible. Indeed, the computational cost of FQ-DCT model as we implemented in LAMMPS code is about three to four orders of magnitude less than the FPMD computational cost for the simulations we performed here.

In order to take into account charge transfer effects explicitly (between water molecules and ions, as well as among water molecules) in classical molecular dynamics simulations, we implemented FQ-DCT model by Rick and co-workers^22–25 into a highly parallel LAMMPS code.^26,35 Some of the parameters for the potentials in systems containing K⁺ were modified compared to the original ones from Rick and co-workers. The new parameters are provided in Table S1 of supplementary material, and Figures S1 and S2 show comparisons.⁵¹ In addition to the charge transfer effect, polarizable effects of water and ions are also taken into account by fluctuating charge¹⁹ and Drude oscillator methods,¹⁷ respectively. In each time step of MD simulations, the amount of charge transfer as a function of relative atomic positions of ions and water molecules is determined first. Then, a self-consistent iteration is performed, so that the potential energy is minimized with respect to the charges on atoms in water molecules and the positions of Drude oscillators attached to ion atoms. We set the convergence criterion for the self-consistency to 1.0 × 10⁻⁵ kcal/(mol Å). The Lennard-Jones and the short-range part of Coulomb interactions are calculated with a cutoff of 10.0 Å. A Particle-Particle Particle-Mesh (PPPM)³⁶ algorithm is used for the long-range Coulomb interaction with the accuracy of 1.0 × 10⁻⁵ kcal/(mol Å). Atoms in water molecule are rigid in these simulations, and SHAKE algorithm³⁷ is used to constrain the bond lengths and angles with tolerance of 1.0 × 10⁻⁴. The NVT ensemble at the temperature of 298 K was generated by Nose-Hoover thermostat³⁸ with a damping parameter of 1000 fs. The number of ions and water molecules as well as simulation cell volumes at different salt concentrations is determined from experimental density.³⁹ In each simulation, a 100 ps equilibration run is performed, followed by a 1 ns of a production run (0.5 fs time step).

For the FPMD simulations,⁴⁰ we followed the same procedure/protocols as in the recent work by Ding et al.²¹ CP2K code⁴¹ was used for these simulations. In these simulation, hybrid Gaussian and plane waves and its extension with augmented plane waves (GAPW) schemes were used.⁴² The GAPW scheme was applied to Na⁺, K⁺, and Cl⁻ ions to obtain well-converged forces. The plane wave cutoff was set to 280 Ry. Revised Perdew-Burke-Ernzerhof functional (revPBE)^43,44 was used for the exchange and XC functional. Grimme D3 type empirical dispersion correction⁴⁵ (revPBE-D3) was included to account for the dispersion interaction. Supplementary material⁵¹ contains a validation of this particular XC approximation against the more appropriate range-separated approximation of LC-BLYP.^46,47 Goedecker-Teter-Hutter pseudopotentials⁴⁸ were used for the core electrons; for Na⁺ and K⁺, the semicore electrons were included as valence electrons. Double-zeta split valence basis sets⁴⁹ were used for all atomic species as in the work by Ding et al.²¹ The simulation was performed with a time step of 0.5 fs, and the Kohn-Sham wavefunctions were optimized at each step with the orbital transformation method⁵⁰ to ensure the system follows the Born Oppenheimer potential energy surface. The NVT ensemble at the temperature of 298 K was generated by Nose-Hoover thermostat.³⁸ Additional computational details and parameters are provided in the supplementary material.⁵¹ 

Y.K. and Y.Y. gratefully acknowledge support by the donors of the Petroleum Research Fund, Grant No. 52494-DNI6 and National Energy Research Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 for computational resources. We thank Professor S. Rick for providing to us the parameters for FQ-DCT model (Ref. 22). We also thank Dr. Yun Ding (ETH, Zurich) for providing us the details of first principles molecular dynamics simulations presented in Ref. 21.