Common Cations Are Not Polarizable: Effects of Dispersion Correction on Hydration Structures from Ab Initio Molecular Dynamics

We employed density functional theory-based ab initio molecular dynamics simulations to examine the hydration structure of several common alkali and alkali earth metal cations. We found that the commonly used atom pairwise dispersion correction scheme D3, which assigns dispersion coefficients based on the neutral form of the atom rather than its actual oxidation state, leads to inaccuracies in the hydration structures of these cations. We evaluated this effect for lithium, sodium, potassium, and calcium and found that the inaccuracies are particularly pronounced for sodium and potassium compared to the experiment. To remedy this issue, we propose disabling the D3 correction specifically for all cation-including pairs, which leads to a much better agreement with experimental data.

In order to demonstrate for cations the artificial effect of the dispersion correction, we first present gas-phase interaction energies (E int ) of cation-water and cation-cation pairs calculated using the revPBE and revPBE-D3 density functionals, as compared to the coupled cluster [CCSD(T)] method as the "golden standard" (Figure 1).At large to intermediate separations, the interaction energies are determined primarily by charge-dipole and charge-charge electrostatic interactions in the case of cation-water and cation-cation pairs, respectively.These interactions are well captured by DFT, hence the very good agreement between the revPBE and CCSD(T) potential energy curves.Note that the revPBE method does not account for dispersion interactions while CCSD(T) does, confirming that dispersion is negligible for the cations investigated here.
The D3 correction incorporated in the revPBE-D3 functional results for the investigated systems in a spurious shoulder on the potential energy curves (Figure 1).This artificial stabilization of 10-30 kJ/mol, which is a direct consequence of the D3 parametrization treating in terms of dispersion atomic cations as neutral atoms, indeed occurs around the positions of hypothetical van der Waals minima between the corresponding neutral pairs.Note that the D4 parameterization, which considers not only the atomic specificity but also the oxidation state, does not remove -but only slightly reduces -the size of this artifact (see SI for more details).Also note that in relative terms, the artifact is larger for sodium and potassium than lithium and calcium (see SI) due to weaker electrostatic interactions in the case of the former cations.Finally, we stress that the purpose of the present gas phase calculations is to point directly to the artifact due to the inappropriate use of the atomic dispersion term for the corresponding cation, which is helpful for understanding the analogous effect on the structural properties of aqueous salt solutions in the condensed phase (vide infra).
Two sets of AIMD simulations of a single aqueous Li + , Na + , K + or Ca 2+ (see SI) cations were performed: one with the D3 dispersion correction enabled for all atoms (i.e., default setting) and second with the D3 term disabled for all atomic pairs involving the cation.The hydration structure of the cations was quantified using radial distribution functions (RDFs) and running coordination numbers (RCNs) of the surrounding water molecules, as presented in Figure 2. We see sizable artifacts due to the D3 correction for sodium and potassium cations.Namely, the presence of the dispersion term on these cations leads to a looser hydration shell which is shifted to a larger separation from the cation and contains about one additional water molecule.The artificially large dispersion term thus effectively increases the interaction of the cation with the surrounding water molecules but at the same time reduces its ability to reorient and thus bring closer the neighboring water molecules.This is because dispersion interactions depend much less on the ion-water molecule relative orientation than the electrostatic ion-dipole interactions.Comparison between the four investigated cations shows that the strength of the electrostatic interaction plays an important role.In particular, for the two cations with a high charge density and thus stronger electrostatic interactions, i.e., for lithium and calcium, the relative importance of the dispersion interaction artifact is smaller than for sodium and potassium, which have a lower charge density (Figure 2).Quantitatively, as detailed in Table I, upon zeroing the cationic D3 terms, a shift of 0.05 Å, 0.10 Å in the position of the first peak of the Na-O and K-O RDF respectively is apparent as well as significant changes in the average coordination number (CN), moving from 6.04 to 5.32 for sodium and 7.70 to 6.39 for potassium.It is reassuring that disabling the cationic D3 correction leads to better agreement with neutron scattering experiments for sodium and potassium (Table I), in particular in terms of the average coordination numbers.Note that accurate hydration structure was also achieved by the SCAN functional, 14 which bypasses the use of empirical dispersion being, however, less accurate for the structure of liquid water compared to revPBE-D3.
Next, we simulated concentrated (4 M) aqueous solutions of lithium, sodium, and potassium chloride with and without the cationic D3 correction which allowed us to acquire cation-water, cation-cation, and cation-anion RDFs.These conditions are close to the experimental setup used in neutron scattering measurements, which are typically conducted at higher salt concentrations, compared to the infinite dilution.The cation-water RDFs shown in the top part of Figure 3  play the same characteristics as those from the infinitely dilute systems, but with the CN of sodium cation closer to the experimental value found in Table I.In concentrated solutions, the second peak of the Na-O RDF is also affected by the D3 correction on sodium, in contrast to the infinitely diluted system.However, the general trend of increased CNs for cations with the D3 correction still holds.
The D3 correction has a pronounced effect on sodium and potassium cations, while there are also non-negligible changes for the lithium cation.Namely, the first peaks of the cationwater RDFs are larger when the D3 correction is applied to the cations.This increase is due to the presence of the shoulders at 3-4 Å at the potential energy curves between cations and water, as well as at 4-5 Å at the potential energy curves between a pair of cations (Figure 1).
Simulations of concentrated solutions provide insights also into cation-cation interactions, as well as interactions between cations and anions, as shown in the middle and bottom panels of Figure 3. Again, the peak intensities are larger when the D3 is applied to cations.This is the most pronounced for sodium, and in the case of K-K RDF, there is also a noticeable shift in the position of the first peak.Similar, although less noticeable, effects are also observed in the Li-Li RDF.
The cation-anion RDFs, shown in the bottom panels of Figure 3, exhibit significant differences with or without dispersion.The first peak is more pronounced when the D3 correction is turned off for metallic cations, and in the case of the Na-Cl and K-Cl RDFs and the position of the first maximum and first minimum are shifted by about 0.2 Å, 0.3 Å respectively.These RDFs can be directly converted to the free energies of ion pairing (Equation 1), as shown in Fig- ure 4. The above observations have important implications for the parametrization of classical ionic force fields, where the AIMD-based free energies of ion binding have been used as a reference. 18,19nally, we compared our simulated results for a 4 M potassium chloride solution to neutron scattering data. 17The total radial distribution function (RDF) for potassium was obtained

K Cl
D3 for all atom types D3 excluded for K FIG. 4. Free energy profiles of cation-chloride pairing as a function of distance for lithium chloride (top, green) and sodium chloride (middle, blue) and potassium chloride (bottom, red) obtained from the corresponding RDFs from Figure 3.
as a weighted sum of simulated RDFs for K-O, K-K, K-H, and K-Cl with appropriate weights from the reference. 17We then obtained its reciprocal space equivalent via Fourier transform in order to directly compare it with experimental neutron scattering data.Our simulated RDFs, with and without the D3 correction for potassium, are shown along with the experimental curves in Figure 5 in both reciprocal and real space.Clearly, better agreement with experiment was obtained without the D3 correction for potassium cation.
In this study, we used DFT-based AIMD simulations to examine the hydration structure of common alkali metal cations.We found that the commonly used D3 correction scheme, which does not take into account the particular oxidation state, leads to inaccuracies in the hydration structure, as it assigns cations with dispersion coefficients appropriate for their neutral counterparts.This results in highly polarizable and large cations which does not reflect experimental reality.We evaluated this effect for cations including lithium, sodium, and potassium both at infinite dilution and for 4 M chloride solutions, which is particularly relevant when compared to the neutron scattering experiments.To fix this issue, we propose to disable the D3 correction specifically for cation-including pairs, while still keeping the D3 term for solvent molecules (and anions).By comparison to neutron scattering data, we show that this simple correction significantly improved the hydration structures of the cations, in particular of sodium and potassium.Therefore, we conclude that the commonly used D3 dispersion correction is often applied incorrectly when soft single-charged metallic cations are present.
As the computational power continues to increase, AIMD simulations of salt solutions are becoming a standard tool of computational chemistry.These simulations offer good accuracy when condensed systems are addressed and may provide benchmark data for less accurate methods such as molecular dynamics with empirical potentials or data-driven methods, e.g., neural network potentials.In this context, it is of key importance to properly account for dispersion interactions when studying alkali metal cations, which can be easily accomplished by zeroing the D3 term for all pairs involving a cation.

COMPUTATIONAL DETAILS
The present simulations of bulk aqueous solutions were performed using the CP2K 9.1 package with the Quickstep module 20 in the canonical ensemble at 300 K employing three-dimensional periodic boundary conditions.Energy and forces were calculated on the fly using the revPBE 21,22 density functional.The dispersion interactions were accounted for using the D3 correction 8 scheme using two-body terms and zero damping and were either enabled or disabled for all pairs involving cations.Using the hybrid Gaussian and plane wave approach, 23 Kohn-Sham orbitals were represented by a TZV2P basis set 24 and the charge density by an auxiliary plane-wave basis up to a 400 Ry cutoff.The core electrons were replaced by Goedecker-Teter-Hutter (GTH) pseudopotentials. 25o simulate each system using AIMD, we first preequilibrated it by means of force-field molecular dynamics (FFMD) in GROMACS 2022.2 26 in the canonical ensemble for 10 ns with the SPC/E water and scaled-charge ion models (Reference 27-30 for Li + , Na + , K + and Ca 2+ respectively).Five initial structures were extracted from the FFMD trajectory with a 2 ns stride for each system which served as independent starting points for the subsequent AIMD simulations.These five structures were then equilibrated for 5 ps using the Langevin thermostat 31 with a 50 fs time constant.Next, the stochastic velocity rescaling (SVR) thermostat 32 was employed with a 1 ps time constant, and 20 and 40 ps of production runs were acquired for each structure.This amounts to a total simulation time of 100 and 200 ps per system for the dilute and 4 M solution with chloride counterions.Simulations with and without the cationic D3 correction were started from the same initial structures.
The dilute systems consisted of one cation surrounded by 128 water molecules in a cubic box.The box size was chosen to match the density of pure liquid water of 1 kg/dm 3 resulting in 15.656 Å for lithium, 15.692 Å for sodium, 15.728 Å for potassium and 15.730 Å for calcium cation.The 4 M solutions contained 8 cations, 8 chloride anions, and 111 water molecules.The box sizes for these systems, as determined based on the experimental densities of the salts, were 15.314 Å for LiCl, 33 15.328 Å for NaCl 34 and 15.524 Å for KCl. 35ree energy profiles of ion pairing (A) were calculated from the RDFs in a standard manner as follows where k B is the Boltzmann constant and T temperature.Finally, gas-phase calculations were performed using the ORCA 5.0.3 software. 36,37The interaction energy (E int ) curve was obtained as an energy scan of the distance between the cation and the oxygen of the water molecule or another cation, as described by the equation: where the subscript denotes the system (M: metal cation, X: water or metal cation), and the superscript indicates the basis set used.Such E int calculation includes a counterpoise correction for the basis set superposition error. 38The revPBE 21,22 density functional (optionally equipped with the D3 8 or D4 11,12 dispersion correction) and coupled cluster single, double and perturbative triple excitations [CCSD(T)] 39 methods were used for the E int calculation.In all the gasphase calculations, the def2-TZVPP 40  In Figure S1, we show gas-phase interaction energies of a cation with water (left panel) and the same cation (right) obtained according to the Equation 2 in the main text for Li + , Na + , K + , and Ca 2+ .E int calculated by the revPBE density functional and the CCSD(T) (in the case of calcium) method are compared on the left-hand side y-axis, while the isolated D3 and D4 dispersion correction to the revPBE energies on the right-hand side y-axis.It is worth a note that we subtracted Coulomb potential from the cation-cation curves on the righthand side panels of Figure S1 in order to facilitate the readout of features that would be otherwise covered by the chargecharge repulsion.

FIG.
S1. E int of cation-water (left) and cation-cation (right) as a function of the distance for lithium, sodium, potassium, and calcium cations from top to bottom.E int obtained at revPBE level is plotted on the left y-axis in red and the D3 (solid) and D4 (dashed) dispersion correction on the right y-axis in grey.Additional CCSD(T) (dashed green line) result is displayed for calcium cation.Note that Coulomb interaction potential was subtracted from the cation-cation curves.a) Electronic mail: hseara@gmail.comb) Electronic mail: pavel.jungwirth@uochb.cas.czThe D3 contributes to the E int at distances larger than where the minimum is located, which agrees with the results presented in Figure 1 of the main text.This contribution amounts roughly to 10 kJ/mol for cation-water curves.For the cationcation cases, it increases to 20-30 kJ/mol.The magnitude of the D4 dispersion correction is generally smaller than that of D3.However, the problematic position of the peak is practically the same for both corrections in spite of D4 reflecting better varying electronic structure of the involved species with respect to D3.In the case of calcium, used here as an example of a divalent cation, we see a significant discrepancy between E int calculated either with revPBE and the CCSD(T).The difference is a consequence of the spurious overdelocalization of electrons in the gas phase, which is not surprising for the GGA DFT.This artifact does not occur at the CCSD(T) level.
Additional RDFs of water and chloride anion are shown in this section for the infinitely diluted and concentrated systems in Figure S3.It reveals that neither structure of water nor the hydration of chloride anion is perturbed by the exclusion of metal cation from the D3 dispersion calculation.

S2. CALCIUM CATION AT INFINITE DILUTION
We simulated a single calcium cation at infinite dilution according to the protocol in the Computational Details of the main text.We extracted the Ca-O RDF and RCN with and without the D3 correction applied to the calcium cation.Both setups yielded similar results, with the first peak at 2.43 Å.This agrees well with the value of 2.38 Å obtained in the neutron scattering experiment.S1 The main difference between the two simulations is the average CN.In the case where the D3 correction was applied to the calcium cation, the CN is slightly larger at 6.68 water molecules in the first hydration shell.In contrast, a CN of 6.33 is obtained when the cationic D3 correction is not used.This observation is consistent with the results presented for lithium, sodium, and potassium in the main text.

FIG. 1 .
FIG.1.Interaction energy curves in the gas phase as a function of distance for cation-water (left) and cation-cation (right).The top panel illustrates the orientation of the water molecule with respect to the cation (blue) and the distance r used in the interaction energy scan.The curves are color-coded to indicate the employed method: grey for revPBE-D3, red for revPBE, and dashed green for CCSD(T).Note that the E int is aligned to zero at the largest distance for all cases.

FIG. 2 .
FIG.2.Radial distribution functions of cation-oxygen (full) and their corresponding running coordination numbers (dashed) including the D3 correction for all atoms (grey); and all atoms except Li (green), Na (blue), and K (red) cations.The panels are arranged in the order of lithium (top), sodium (middle), and potassium (bottom).

FIG. 3 .
FIG.3.Radial distribution functions of metal cations in 4 M LiCl (left), 4 M NaCl (middle), and 4 M KCl (right) solutions.The top, middle, and bottom panels show the radial distribution functions of water oxygens, metal cations, and chloride anions, respectively, surrounding metal cations.The results are shown with the D3 correction turned on for all atom types in grey and for all types except Li (green), Na (blue), and K (red) cation.

FIG. 5 .
FIG. 5. Simulated structure factor S (top) and radial distribution function (bottom) for potassium with water and chloride, shown in red (without D3 correction on potassium) and grey (with D3 correction on potassium).The experimental neutron scattering structure factor (top) and the corresponding RDF (bottom) from Reference 17 are depicted in orange for comparison.
FIG. S2.RDFs of calcium-oxygen (full) and corresponding RCN (dashed) including the D3 correction for all atom (grey) and for all kinds but calcium (pink).
FIG. S3.RDFs in the infinitely diluted system (left panels) and in the 4M chloride solution (right panels).Oxygen-oxygen RDFs are displayed for all systems, and oxygen-chloride RDFs only for the 4 M solutions.Cations are color-coded as green (lithium), blue (sodium), and red (potassium).

TABLE I .
Metal-water oxygen radial distribution function first peak positions (RDF max ) and average coordination numbers (CN) for systems at infinite dilution (Figure2) and at 4 M chloride solution (Figure3).Experimental values were adopted from Reference 15 for lithium, 16 for sodium and, 17 for potassium.The use of the cationic D3 correction is indicated in parentheses.
basis set was employed.Supporting information for: Common Cations are not Polarizable: Effects of Dispersion Correction on Hydration Structures from Ab Initio Molecular Dynamics Vojtech Kostal, Philip E. Mason, Hector Martinez-Seara*, a) and Pavel Jungwirth* b) Institute of Organic Chemistry and Biochemistry of the Czech Academy of Sciences, Flemingovo nám.2, 166 10 Prague 6,