Structure and dynamics of aqueous solutions from PBE-based first-principles molecular dynamics simulations

Density functional theory (DFT)-based first-principles molecular dynamics (FPMD) simulation^1–3 has been a powerful approach for accessing physical and chemical processes in biology, chemistry, and materials science, such as electrocatalytic water splitting for solar-to-fuel conversion.^4–9 Despite this success, establishing an accurate and predictive FPMD framework for the simulation of aqueous solutions remains a great challenge. A major difficulty lies in the limitation of the generalized gradient approximation (GGA) that is typically utilized for the exchange-correlation functional in the application of FPMD simulations. In particular, it is well-known that GGA approximations suffer from the self-interaction error¹⁰ and the lack of van der Walls (vdW) interactions,^11–16 which play important roles in the description of water hydrogen bonds. In addition, most FPMD simulations neglect nuclear quantum effects (NQEs) that are critical to account for the quantum nature of light hydrogen atoms in liquid water.^17–20 The interplay between limitations of the GGA approximation and the lack of NQEs is known to lead to several discrepancies between simulations and experiments, including the overstructuring of liquid water commonly found in GGA-based FPMD simulations.²¹ 

A number of approaches have been proposed to improve the description of aqueous solutions in DFT-based FPMD simulations. For example, the inclusion of vdW corrections is known to improve certain properties such as the density of liquid water.²² However, the extent to which these corrections affect the structure and dynamics of aqueous solutions remains unclear, as it often depends on the approximation employed for the vdW interaction.^22–28 In addition, hybrid exchange-correlation functionals can be utilized to alleviate the self-interaction error in GGA functionals,²⁹ and several approaches have been proposed to include NQEs in the study of aqueous solutions.²⁰ Nevertheless, these simulations are computationally intensive, and thereby limited to short time scales and relatively small system sizes.^20,30–34 A high computational cost has also prohibited the routine application of these advanced theoretical methodologies to a wide variety of systems, including aqueous interfaces.

Other approaches for improving the description of aqueous solutions include the use of sophisticated force fields based on high level quantum chemistry methods and DFT.^35–37 For example, recent work by Paesani and co-workers has demonstrated that force fields derived from many-body expansion of the interaction energy are able to effectively reproduce static and dynamical properties of gas and liquid phases of water.^35,38,39 While these developments represent an essential step in establishing a predictive framework to study aqueous solutions, their applications to solutions involving heterogeneous components remain a significant challenge.⁴⁰ 

On the other hand, given the applicability to a wide range of systems and a good balance between the accuracy and computational expense, GGA-based simulations remain an active field in the study of aqueous solutions. These simulations often rely on the observation that the GGA inaccuracies in the description of water structure and hydrogen bonds at ambient conditions can be artificially “corrected” by performing simulations at elevated temperatures.^41–46 For example, it is known that a temperature around T = 400 K can be used for the Perdew-Burke-Ernzerhof (PBE) functional⁴⁷ to qualitatively recover water structure at ambient conditions.^41,42 This computational framework has also been extended to investigate complex aqueous systems, such as those with solvated ions;^{31–33,48–50} however, a detailed understanding of how the choice of elevated temperatures influences properties of different solution components is largely lacking. In particular, it remains unclear how the simulation temperature affects the description of the ion solvation structure in salt solutions, and if the use of an elevated temperature would lead to a good description of the dynamical properties of these solutions, as has been shown for pure water.

The aim of this paper is to provide a systematic study of temperature effects on the structure and dynamics of liquid water and complex aqueous solutions in GGA-based simulations. We chose to study three aqueous solutions with representative ions, including Na⁺, K⁺, and Cl⁻. We focused on the PBE exchange-correlation functional, and we performed simulations at seven different temperatures in the range of 300–460 K for each solution. For liquid water, we show that the use of elevated temperatures between 390 and 400 K enables a good description of water structure and diffusivity, consistent with previous studies.^41,42 Furthermore, we show that this computational setup is transferable to complex aqueous solutions, providing ion solvation structures and ion effects on water dynamics consistent with experiments. Overall, our results suggest that an elevated temperature around 390–400 K with the PBE functional can be used for the description of structural and dynamical properties of liquid water and complex solutions with solvated ions at ambient conditions.

The rest of the paper is organized as follows. We describe our computational methods in Section II. We re-examine the temperature effect in the simulation of liquid water, focusing on the oxygen-oxygen radial distribution function (RDF), self-diffusion coefficient, and water dipole moment, all of which are available from experiment in Section III. We also discuss the time scale effects on the structural and dynamical properties of liquid water. We then investigate the temperature effects in the simulation of solutions with solvated ions in Section IV, focusing on the ion solvation structure and specific ion effects on water dynamics. Concluding remarks are given in Section V.

Aqueous solutions with solvated ions were modeled by cubic cells consisting of one ion and 63 water molecules, with the cell size chosen to match the experimental density of liquid water at ambient conditions. Counter-ions were not considered in our simulations, and a uniform background charge was employed to neutralize the systems. Simulations of liquid water were carried out using 64 water molecules with identical supercell size. The choice of system size was justified based on classical molecular dynamics simulations of models consisting of one ion and a varying number of water molecules (63, 127, 255, and 511). In particular, we carried out classical molecular dynamics simulations with the CHARMM package,⁵¹ using simple point charge models for water molecules (SPC/E)⁵² and ions.⁵³ We found that the ion-oxygen radial distribution functions (RDFs) of these simulations are not sensitive to the system size and indicate that a system containing 63 water molecules is sufficient for a proper description of the first and second solvation shells of the ions.

Born-Oppenheimer MD (BOMD) simulations were carried out using the Qbox code⁵⁴ with interatomic forces derived from density functional theory (DFT) and the PBE exchange-correlation functional. The Brillouin zone is sampled at the Γ–point. The interaction between valence electrons and ionic cores was represented by norm-conserving pseudopotentials, and the electronic wave functions were expanded in a plane-wave basis set truncated at a cutoff energy of 85 Ry. We used Hamann-Schluter-Chiang-Vanderbilt (HSCV) pseudopotentials⁵⁵ for hydrogen, oxygen, and chlorine atoms. Troullier-Martins and optimized norm-conserving Vanderbilt pseudopotentials^56,57 were employed for sodium and potassium, respectively. The 2p semicore state of sodium and the 3s and 3p of potassium were treated as valence states. All hydrogen atoms were replaced with deuterium in order to maximize the allowable time step, which was chosen to be 10 a.u. in all simulations.

To examine temperature effects on the aqueous solutions, we carried out simulations at constant temperatures (NVT conditions) ranging from 300 K to 460 K (300, 340, 380, 400, 420, 440, and 460 K) for all aqueous solutions, using a velocity scaling thermostat. For each solution, we collected statistics over 50 ps BOMD simulations after 20 ps equilibration runs.

In order to test the sensitivity of our results to time scale effects, we extended the NVT simulation at 380 K up to 170 ps and we collected statistics for 150 ps after 20 ps equilibration. Fig. 1 shows the g_OO(r) computed for different simulation lengths (10 ps, 30 ps, 50 ps, and 150 ps), indicating only a minor difference in the peak intensity between 30 ps and 150 ps trajectories. In addition, the g_OO(r) intensity computed for the 50 ps simulation coincides with that obtained with the 150 ps trajectory within 0.05, indicating that g_OO(r) results are converged for the simulation length of 50 ps. On the other hand, g_OO(r) obtained with only 10 ps exhibits significant errors up to 0.2, particularly around the first maximum.

Similar to g_OO(r), the water diffusion coefficient (D_w) also converges for the simulation length of 50 ps. Specifically, we computed D_w from the slope of the oxygen mean square displacement (MSD) in the range of 1–10 ps, which was obtained by averaging over the trajectories with multiple starting configurations evenly spaced by 0.12 ps. We obtained a value of 1.5 and 1.3 (×10⁻⁵ cm²/s) for D_w using a trajectory of 30 ps and 50 ps, respectively, in agreement with the value of 1.4 × 10⁻⁵ cm²/s computed for the full simulation within a few percent. Therefore, we conclude that 50 ps is sufficient to converge structure and diffusivity of liquid water with respect to the simulation time.

We first re-examine temperature effects on the water structure with the focus on the oxygen-oxygen radial distribution function (g_OO(r)). As expected and demonstrated in Fig. 2, water structure shows a strong dependence on the simulated temperature. In addition, Table I presents details of the g_OO(r) computed at different simulation temperatures, including the intensity and peak position of the first maximum (g_max,1, r_max,1), the first minimum (g_min,1, r_min,1), and the second maximum (g_max,2, r_max,2). We find that the simulation at 300 K leads to overstructured liquid water, with the intensity of the first peak overestimated by more than 40.0%. On the other hand, the use of an elevated temperature around 400 K yields peak intensities in qualitative agreement with experiments (Fig. 3, top panel). Specifically, the average mean absolute error (MAE) of peak intensities (g_max,1, g_max,2, and g_min,1) computed for simulations at 400 K is 10.5% and 8.5% compared to experimental data reported in Refs. 58 and 59, respectively.

In contrast to the intensity, g_OO(r) peak positions are not sensitive to the simulation temperature, showing good agreement with experiment in a wide temperature range of 300–460 K (Fig. 3, lower panel). Specifically, the MAEs computed for the peak positions (r_max,1, r_max,2, and r_min,1) vary within 3.6%–1.2%, with the errors of less than 0.2 Å for all simulation temperatures. In addition, the first-neighbor coordination number (n_OO) behaves in the same manner, showing a weak dependence on the simulation temperature, and yields a value of 4.14 at 400 K. Although experimental data of n_OO exhibit a large uncertainty,^58,59 the theoretical result obtained at 400 K, where the peak intensities are in good agreement with experiment, is consistent with a more recent experimental measurement at ambient conditions within 0.15.

Collectively, we conclude that the use of elevated temperatures around 400 K generally provides a good description of ambient water structure. It is also worth pointing out that the intensity of the first and second maxima decreases linearly with temperature, but with a slightly smaller slope above 420 K (Fig. 3, top panel), indicating a subtle change in the water structure around 420 K associated with the collapse of the second solvation shell and a rapid increase in the first shell coordination number n_OO above 420 K, as shown in Table I. Finally, we note that adjusting the simulation temperature may not be sufficient to recover the experimental g_OH(r) which requires the explicit inclusion of NQEs.³⁴ 

Similar to structural properties, water dynamics also exhibits significant dependence on the simulation temperature. This is demonstrated in Fig. 4 that shows a quadratic dependence of the water diffusion coefficient on the simulation temperature. When compared to experimental values at ambient temperature, the diffusion coefficient D_w is underestimated by a factor of 13 for the simulation at 300 K, whereas it is overestimated by at least a factor of 1.4 for those performed at T ≥ 420 K. Our calculations show that the best agreement with experimental data is obtained for the simulation performed at 400 K.

In order to make a quantitative comparison between the computed diffusion coefficients and experiments, it is important to take into account the finite size effects that are present in the simulations.^64–67 We estimated this correction using the expression proposed by Dünweg and Kremer,^64,65

where k_B is the Boltzmann constant, T is the temperature, ζ is a parameter equal to 2.837 29, η is the viscosity of liquid water whose value was taken from experimental data,⁶⁸ and L is the size of the cubic simulation cell. When ΔD_w is taken into account, we find that PBE simulations at 390 K are in close agreement with the experimental value of D_w at ambient conditions. Accordingly, our results indicate that PBE simulations performed in the temperature window of 390–400 K yield water diffusivity in good agreement with experiments. We note that the use of theoretical water viscosity is likely to lead to a smaller correction due to the underestimation of the water viscosity in PBE simulations,⁶⁶ and therefore does not affect our conclusions.

Our calculations also indicate that liquid water obtained with PBE simulations below 340 K is in a glassy state with a low diffusivity; particularly, while D_w computed at 400 K is comparable to experimental data, D_w obtained at 340 K is about eightfold smaller. This observation is consistent with the previous studies of PBE water at this temperature range.^45,66 It is worth mentioning that calculations of D_w in this state are known to suffer from slow convergence and are quite sensitive to the initial conditions.⁴⁴ Finally, we note that our discussion is based on calculated diffusion coefficients; on the other hand, the PBE melting temperature of ice was estimated to be 417 K,⁶⁹ indicating that PBE water at 400 K is likely to correspond to a supercooled state.

We investigated temperature effects on the electronic properties of liquid water, focusing on the water dipole moment and band gap of liquid water. We find that these quantities are also affected by the simulation temperature, although the effects are less significant compared to the structural and dynamical properties. In particular, the distribution of the molecular dipole moment (μ) obtained with maximally localized Wannier functions (MLWFs) shows a shift toward smaller values (Fig. 5), yielding an average reduction of 0.3 D (9.2%) when the temperature increases from 300 K to 460 K. Although the molecular dipole moments obtained at different temperature are within the error bar of the experimental value of 2.9 ± 0.6 D,⁶⁰ the use of an elevated temperature generally leads to better agreement with the mean experimental result, consistent with the trend observed for structural and dynamical properties. In particular, the water dipole moment yields a value of 3.04 D at 400 K, in agreement with the average experimental result within 5%.

In addition, we investigated the temperature effect on water band gap and band edges, which play an important role in the determination of the efficiency of solar-to-hydrogen conversion devices.⁹ We find that the DFT band gap weakly depends on the temperature, showing an average reduction of 0.4 eV when the temperature is raised from 300 K to 400 K. This band gap closing is a result of the variation in both water valence band maximum (VBM) and conduction band minimum (CBM), as demonstrated in Fig. 6 for 100 configurations extracted equally in time from 300 K to 400 K trajectories. Specifically, comparison between these two trajectories shows that the VBM and CBM exhibit average variations of 0.2 eV. Such variations in the band gap and band edge positions reflect the interplay between structural and electronic properties and may influence the results obtained with high level theories, such as the G₀W₀ approach or DFT with more sophisticated approximation for the exchange-correlation functionals.^70–73 We emphasize that direct comparison with experiment for the water band gap and band edges is beyond the scope of this work, as it requires the use of higher level theories beyond DFT/PBE.^74–77 

To summarize, our results are consistent with previous studies,^41,42,45,66 showing that simulations with the PBE functional at temperatures 390–400 K allow for good description of structure and dynamics of liquid water at ambient conditions. It is important to notice that even when an elevated temperature is used, discrepancies with experiments may remain for certain properties due to fundamental limitations of the PBE exchange correlation functional and the neglect of NQEs. For example, it is well-known that a sizable red shift of about 200 cm⁻¹ relative to experimental data remains in water IR vibrational spectra even for simulations carried out at 400 K.^78–80 

We now turn to discuss solvated ions in liquid water. We first investigate the ion solvation structure, focusing on the ion-oxygen radial distribution function and the comparison with corresponding experimental results. In this case, the position of the first maximum (r_XO,1) in the ion-oxygen radial distribution function (g_XO(r)) and the oxygen coordination number (n_XO) in the first ion solvation shell are often utilized for comparison with experiments, due to the lack of available experimental data for the g_XO(r) peak intensity.

Similar to liquid water, our simulations show that increasing the simulation temperature also results in a reduction in the ion-oxygen RDF peak heights (Fig. 7). In addition, we report in Figs. 8 and 9 direct comparison between experiments and theoretical results for r_XO,1 and n_XO, respectively. We find that r_XO,1 exhibits a weak dependence on the simulation temperature, showing results in good agreement with experimental values of 2.34–2.43 Å,^81,82 2.65–2.73 Å,^81,83 and 3.11–3.16 Å^81,83 for Na⁺, K⁺, and Cl⁻, respectively. Similarly, the computed n_XO is in agreement with experimental data for the whole range of the simulation temperature. Such an agreement partly stems from a large variation in experimental data for the coordination number, i.e., 5.3 ± 0.8,⁸¹ 6.0–6.1 ± 1.0,^81,83 and 6.4–6.9 ± 1.0^81,83 for Na⁺, K⁺, and Cl⁻, respectively. We note that the computed oxygen coordination of Cl⁻ at 300 K is slightly smaller than the lower bound of the experimental data; however, this discrepancy between theory and experiment is smaller than the standard deviation of theoretical calculations, which is 0.44 as obtained from averaging over 10 segments of 5 ps long of the Cl⁻ solution simulation. Overall, our results indicate that agreement with experimental data of r_XO,1 and n_XO could be achieved with any simulation temperature in the range of 300-460 K. We therefore conclude that a simulation temperature of around 390-400 K can be used to simultaneously provide structure of water and solvated ion consistent with experiments.

Although the choice of temperature has only a small effect on n_XO and r_XO,1, it may significantly alter the local structure of the ion solvation shell. Specifically, Fig. 9 shows that the coordination number n_XO is generally larger at higher temperatures, particularly above 420 K. The magnitude of this variation reflects the hydration strength of the ions, for example, n_Na−O shows a weaker variation compared to K⁺ and Cl⁻ when the temperature increases from 300 K to 460 K, indicating a stronger and more rigid solvation shell for Na⁺. Accordingly, we conclude that temperature effects on the ion solvation structure vary for different ions, depending on the strength of ion-water interactions; however, such effects are not significant enough to be verified based on the comparison between theory and experiments on n_XO and r_XO,1.

The effect of simulation temperature on ion local structure is also reflected in the distribution of the oxygen coordination number in the ion first solvation shell. In particular, the distribution is broader at higher temperatures (Fig. 10) for all ions, indicating more flexible and disordered ion solvation shells. Evidence of the variation in the local structure of ion solvation shells is also found in broader distributions of the tilt angle between the ion–oxygen distance and the water dipole moment with increasing temperatures (Fig. 11). However, the average tilt angle only shows slight variations within 2.0° in the range of 300-400 K for all ions, again indicating weak temperature effects on the average ion solvation structure.

Finally, we investigate temperature effects on water dynamics in the presence of solvated ions, which are further complicated due to their subtle interplay with ion hydration behavior. In particular, it is known that water dynamics in salt solutions deviates from pure liquid water, and the concepts of structure-breaking and structure-making have been used in order to provide a qualitative explanation of ion effects on water dynamics, such as the ionic specificity in the Jones–Dole viscosity B coefficient.⁸⁴ Specifically, ions with a positive B reflect a strong hydration and are deemed to be “structure making,” and the others as “structure breaking.” In this regard, K⁺ and Cl⁻ are classified as structure breakers, while Na⁺ is a structure maker.⁸⁴ Consistent with this argument, it has been shown experimentally that the water diffusion coefficient in aqueous solutions decreases with increasing NaCl concentration, while it increases with increasing KCl concentration.^85,86

On the other hand, FPMD studies of ion effects on water diffusivity are not conclusive. Simulations using the Hamprecht-Cohen-Tozer-Handy functional⁸⁷ at 300 K show that water diffusivity is enhanced by both Na⁺ and K⁺,⁸⁸ while it was reported to be suppressed by the same ions in a recent FPMD simulation using similar setup.⁸⁹ In addition, ion effects may depend on the choice of the exchange-correlation functional, i.e., simulations with the revPBE functional predict that K⁺ enhances the water diffusivity, while it suppresses this motion when vdW corrections are taken into account.⁹⁰ Likewise, providing an accurate description of ion effects on water dynamics is also a challenge for classical MD simulations. For example, recent studies show that MD simulations with both fixed charge and polarizable force fields failed to qualitatively predict the experimental trends of water diffusivity as a function of the salt concentration.^91,92

To investigate the temperature effect on water diffusivity in the aqueous solutions, we computed the ratio (⁠$Dwa/Dw$⁠) between water diffusion coefficients of the solutions (⁠$Dwa$⁠) and that of pure liquid water (D_w). When this ratio is larger than one, water diffusion in the presence of solvated ions is faster than bulk water, and vice versa. We focused on the simulation temperature window of 380–420 K, at which the computed water diffusion coefficient is in better agreement with experiment. Calculations of water diffusion coefficients in the aqueous solutions were carried out using the same procedure as described in Section III for liquid water.

Our results presented in Table II generally indicate that water diffusivity of PBE aqueous solutions obtained with 380–420 K simulations is enhanced by K⁺ and Cl⁻, whereas it is suppressed by Na⁺. In addition, Table II shows that the magnitude of specific ion effects on water diffusivity is significantly influenced by the simulation temperature, indicating a subtle interplay between ion hydration behavior and temperature effects. In particular, simulations at 380 K exhibit clear ion effects on the water diffusivity, while these effects are less obvious at higher temperatures, particularly at T = 420 K for K⁺ and Cl⁻, and at T ≥ 400 K for Na⁺. Overall, the trend observed in our simulations supports the classification of “structure breaking” and “structure making” of the ions and suggests that the PBE simulations at temperatures around 390–400 K provide a qualitative description of ion effects on water diffusivity. We note that, quantitative estimations of ion effects on diffusion, however, require multiple independent simulations that are beyond the scope of this paper.⁹³ Our results also indicate that the ion effect on water dynamics should be interpreted with care, as it strongly depends on the simulation details, especially the equilibrated temperature.

In conclusion, we present a systematic study of temperature effects in PBE-based FPMD simulations of liquid water and three aqueous solutions with solvated Na⁺, K⁺, and Cl⁻ ions. Consistent with previous studies, we show that a simulation temperature of 390–400 K leads to water structure and dynamics in qualitative agreement with experiments. Furthermore, we demonstrate that this computational setup is transferable to more complex aqueous solutions that include solvated ions. In particular, we find that this computational protocol provides ion solvation structures consistent with experiments, as well as a qualitative description of specific ion effects on water dynamics.

We stress that this work has mainly focused on the structure and self-diffusion of aqueous solutions, and adjusting the simulation temperature may not be sufficient to properly describe other properties, such as IR vibrational spectra and optical properties. In addition, the results presented here were obtained with NVT simulations at the experimental density of liquid water; simulations at an elevated temperature under NPT conditions are likely to lead to poor agreement with experiments for the water density, which is known to be underestimated by the PBE functional under ambient pressure and temperature conditions.²² Finally, several questions remain open in the PBE-based simulations that are the topic of future studies, including temperature effects on ion pair formation.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344; part of this work was supported by DOE/BES (Grant No. DE-SC0008938). T.A.P. acknowledges support from the Lawrence Fellowship. Computational support was from the LLNL Grand Challenge Program. We thank Alex Gaiduk, Francois Gygi, and Giulia Galli for fruitful discussions.