Machine learning surrogate models for particle insertions and element substitutions

The real potential of species X is defined as the free energy difference between X in the gas phase just outside the water surface and X in the water. The real potential is the negative value of the work function of species X and differs from the solvation free energy by the product of the water surface potential and the charge of species X. As in the traditional definition, we define the concentration at standard states as 0.1 MPa (1/24.46 mol L^-1) for the gas phase and 1 mol L^-1 for the liquid phase. As explained later, actual TI computations were conducted using the same periodic cell (with a concentration of 0.87 mol L^-1) for both gaseous and aqueous species. The equations presented in this section are also derived for this condition. After the computations, the change in free energy caused by the difference in concentration was corrected using the ideal gas model.

In the present article, I propose two TI schemes shown in Fig. 1: particle insertion and element substitution. This subsection explains the particle insertion method. The real potential is calculated as an integral of the expectation value of the derivative of the Hamiltonian with respect to the coupling parameter $\lambda$, which specifies a path from the reference state ($\lambda$ = 0) to the interacting system of interest ($\lambda$=1): Kirkwood_JCP_1935 ; Zwanzig_JCP_1954

$\displaystyle\alpha$

$\displaystyle=A_{1}-A_{0}=\int_{0}^{1}\left<\frac{\partial H^{\mathrm{in}}}{\partial\lambda}\right>_{\lambda}d\lambda.$

(1)

The reference state is set to be the non-interacting species X and the solvent without the species X. Along the thermodynamic pathway, the interaction between X and the solvent is gradually turned on. This transition can be realized using the Hamiltonian defined as follows: Blumberger_JCP_2006 ; Dorner_PRL_2018

$\displaystyle H^{\mathrm{in}}$

$\displaystyle=\sum\limits_{i=1}^{N_{\mathbf{a}}}\frac{|\mathbf{p}_{i}|^{2}}{2m_{i}}+\lambda\left(U_{\mathrm{1}}+z\Delta\phi\right)+(1-\lambda)U_{\mathrm{0}},$

(2)

where $\lambda$ is a coupling constant, which specifies a path from the species in vacuum and liquid water ($\lambda=0$) to the mixture of them ($\lambda=1$), $U_{\mathrm{0}}$ and $U_{\mathrm{1}}$ denote potential energies of these two endpoints, $z$ denotes the ionic charge. In any calculation using periodic boundary conditions, the energy changes due to the removal or addition of charged species are not entirely well-defined. To correct for this, the potential gap $\Delta\phi$ is introduced. This potential gap essentially specifies the potential of the vacuum level just outside a water surface, which serves as the common reference point in electrochemistry. The potential gap is determined through a separate slab calculation involving an interface between water and the vacuum.  Cheng_PCCP_2012 ; Ambrosio_JCP_2015 ; Ambrosio_JPCL_2018 ; Jinnouchi_npjComputMater_2024 ; Jinnouchi_arXiv_2024

As done in previous studies, Jinnouchi_PRB_2020 ; Jinnouchi_JCP_2021 ; Jinnouchi_npjComputMater_2024 ; Jinnouchi_arXiv_2024 I decompose the TI equation [Eq.(1)] into two integrals:

	$\displaystyle\alpha$	$\displaystyle=\Delta A_{\mathrm{I}}+\Delta A_{\mathrm{II}},$		(3)
	$\displaystyle\Delta A_{\mathrm{I}}$	$\displaystyle=\int_{0}^{1}\left<\frac{\partial H^{\mathrm{in}}_{\mathrm{I}}}{\partial\lambda_{\mathrm{I}}}\right>_{\lambda_{\mathrm{I}}}d\lambda_{\mathrm{I}},$		(4)
	$\displaystyle\Delta A_{\mathrm{II}}$	$\displaystyle=\int_{0}^{1}\left<\frac{\partial H^{\mathrm{in}}_{\mathrm{II}}}{\partial\lambda_{\mathrm{II}}}\right>_{\lambda_{\mathrm{II}}}d\lambda_{\mathrm{II}}.$		(5)

In the first Hamiltonian $H_{\mathrm{I}}$ [Eq. (4)], as in the previous studies, Timothy_CS_2017 a Gaussian soft repulsive potential $U_{\mathrm{model}}$ is gradually added to the potential energy of the solvent without species X in order to create a space for species X:

$\displaystyle H^{\mathrm{in}}_{\mathrm{I}}$

$\displaystyle=\sum\limits_{i=1}^{N_{\mathrm{a}}}\frac{\left|\mathbf{p}_{i}\right|^{2}}{2m_{i}}+\lambda_{\mathrm{I}}U_{\mathrm{model}}+\sum\limits_{i\notin\mathrm{X}}U_{i}\left(0\right).$

(6)

Here, $N_{\mathrm{a}}$ is the number of atoms in the system, and $U_{i}\left(0\right)$ is the atomic potential energy represented by a MLFF of the solvent without species X. In the second Hamiltonian, the solute-solvent interaction potential $U_{\mathrm{model}}$ is gradually replaced by the interaction potential described by an MLFF:

	$\displaystyle H^{\mathrm{in}}_{\mathrm{II}}$	$\displaystyle=\sum\limits_{i=1}^{N_{\mathrm{a}}}\frac{\left|\mathbf{p}_{i}\right|^{2}}{2m_{i}}+\lambda_{\mathrm{II}}\left[\sum\limits_{i=1}^{N_{\mathrm{a}}}U_{i}(1)+z\Delta\phi\right]+$
		$\displaystyle\left(1-\lambda_{\mathrm{II}}\right)\left[U_{\mathrm{model}}+\sum\limits_{i\notin\mathrm{X}}U_{i}\left(0\right)\right],$		(7)

where $U_{i}(1)$ denotes the atomic potential energy of the solution system involving the species X.

As in previous studies, Bartok_PRL_2010 ; Jinnouchi_PRL_2019 ; Jinnouchi_PRB_2019 the atomic potential energies [$U_{i}\left(0\right)$ and $U_{i}\left(1\right)$] are represented as a linear combination of functions $K\left(\mathbf{x}_{i}\left(\kappa\right),\mathbf{x}_{i_{\mathrm{B}}}\right)$:

$\displaystyle U_{i}\left(\kappa\right)$

$\displaystyle=\sum\limits_{i_{\mathrm{B}}=1}^{N_{\mathrm{B}}}w_{i_{\mathrm{B}}}K\left(\mathbf{x}_{i}\left(\kappa\right),\mathbf{x}_{i_{\mathrm{B}}}\right).$

(8)

Here, the descriptors $\mathbf{x}_{i}\left(\kappa\right)=\mathbf{x}_{i}\left[\rho_{i}\left(\mathbf{r},\kappa\right)\right]$ are the functional of the density distribution function around atom $i$,

	$\displaystyle\rho_{i}\left(\mathbf{r},\kappa\right)$	$\displaystyle=\sum\limits_{j\notin\mathrm{X}}f_{\mathrm{cut}}\left(\left|\mathbf{r}_{j}-\mathbf{r}_{i}\right|\right)g\left(\mathbf{r}-\left(\mathbf{r}_{j}-\mathbf{r}_{i}\right)\right)+$
		$\displaystyle=\kappa\sum\limits_{j\in\mathrm{X}}f_{\mathrm{cut}}\left(\left|\mathbf{r}_{j}-\mathbf{r}_{i}\right|\right)g\left(\mathbf{r}-\left(\mathbf{r}_{j}-\mathbf{r}_{i}\right)\right),$		(9)

where the function $g$ is a smoothed $\delta$-function, and $f_{\mathrm{cut}}$ is a cutoff function that smoothly eliminates the contribution from atoms outside a given cutoff radius $R_{\mathrm{cut}}$. Hence, $U_{i}\left(0\right)$ and $U_{i}\left(1\right)$ represent the local energies without and with the group X, respectively. Using these notations, $U_{0}$ and $U_{1}$ can be written as $\sum\limits_{i\notin\mathrm{X}}U_{i}\left(0\right)$ and $\sum\limits_{i=1}^{N_{\mathrm{a}}}U_{i}\left(1\right)$. It should be noted that the MLFF of the fully interacting system $U_{1}$, trained on the FP data of the periodic system, is also shifted by $-z\Delta\phi$. This shift is corrected in Eq. (7).

The ML-aided particle insertion allows for accurate computation of the real potential. However, it requires many numerical integration grids because the integrand changes steeply along the coupling path from the non-interacting state to the fully interacting state. It should also be noted that $U_{\mathrm{model}}$ needs to be designed so that the system remains within the phase space where the MLFF can interpolate. If the real potential $\alpha^{\mathrm{X}}$ of species X is already given by the atom insertion method, the free energy $A_{1}^{\mathrm{Y}}$ of species Y, which is generated by substituting elements in X with other elements, can be obtained by computing the free energy change $\Delta A_{1}^{\mathrm{X\rightarrow Y}}$ associated with the substitution of the elements in the fully interacting system as follows:

$\displaystyle A_{1}^{\mathrm{Y}}$

$\displaystyle=A_{0}^{\mathrm{X}}+\alpha^{\mathrm{X}}+\Delta A_{1}^{\mathrm{X\rightarrow Y}},$

(10)

where $A_{0}^{\mathrm{X}}$ is the free energy of X in vacuum, which can be easily obtained by the ideal gas model. The free energy change $\Delta A_{1}^{\mathrm{X\rightarrow Y}}$ can be calculated as:

	$\displaystyle\Delta A_{1}^{\mathrm{X\rightarrow Y}}$	$\displaystyle=\int_{0}^{1}\left<\frac{\partial H^{\mathrm{sb}}}{\partial\xi}\right>_{\xi}d\xi.$		(11)
	$\displaystyle H^{\mathrm{sb}}$	$\displaystyle=\sum\limits_{i=1}^{N_{\mathbf{a}}}\frac{|\mathbf{p}_{i}|^{2}}{2m_{i}}+\xi\left(U_{\mathrm{1}}^{\mathrm{Y}}+z^{\mathrm{Y}}\Delta\phi\right)+(1-\xi)\left(U_{\mathrm{1}}^{\mathrm{X}}+z^{\mathrm{X}}\Delta\phi\right).$		(12)

Here, the TI changes the potential from $U_{\mathrm{1}}^{\mathrm{X}}+z^{\mathrm{X}}\Delta\phi$ to $U_{\mathrm{1}}^{\mathrm{Y}}+z^{\mathrm{Y}}\Delta\phi$ while keeping the mass of the species fixed to that of species X. This computation is particularly feasible when the chemical characteristics of species Y are similar to those of species X. The real potential of species Y can be obtained as follows:

$\displaystyle\alpha^{\mathrm{Y}}$

$\displaystyle=A_{1}^{\mathrm{Y}}-A_{0}^{\mathrm{\prime Y}},$

(13)

where $A_{0}^{\mathrm{\prime Y}}$ is the free energy of the species Y in vacuum with the mass of the species X. The value can be also easily obtained by the ideal gas model.

Although both the ML-aided particle insertion and element substitution methods provide the free energy change, they may introduce non-negligible errors due to inaccuracies in the MLFF models. These errors can be corrected through a TI calculation that transitions from the MLFF to the potential obtained by the FP method: Jinnouchi_npjComputMater_2024 ; Jinnouchi_arXiv_2024

	$\displaystyle\Delta A_{\kappa}^{\mathrm{FP-ML}}$	$\displaystyle=\int_{0}^{1}\left<\frac{\partial H_{\kappa}^{\mathrm{FP-ML}}}{\partial\eta}\right>_{\eta}d\eta,$		(14)
	$\displaystyle H_{\kappa}^{\mathrm{FP-ML}}$	$\displaystyle=\sum\limits_{i=1}^{N_{\mathrm{a}}}\frac{\left|\mathbf{p}_{i}\right|^{2}}{2m_{i}}+\eta U_{\kappa}^{\mathrm{FP}}+\left(1-\eta\right)U_{\kappa}^{\mathrm{ML}},$		(15)

where $U_{\kappa}^{\mathrm{ML}}$ and $U_{\kappa}^{\mathrm{FP}}$ represent the potential energies for the state specified by $\kappa$ ($\kappa$=0 and $\kappa$=1 denote the non-interacting and fully interacting states, respectively). These energies are calculated using the MLFF and FP methods, respectively. Utilizing $\Delta A_{0}^{\mathrm{FP-ML}}$ and $\Delta A_{1}^{\mathrm{FP-ML}}$, the free energy change based on the FP method $\Delta A^{\mathrm{FP}}$ is calculated as:

$\displaystyle\Delta A^{\mathrm{FP}}$

$\displaystyle=\Delta A^{\mathrm{ML}}+\Delta A_{1}^{\mathrm{FP-ML}}-\Delta A_{0}^{\mathrm{FP-ML}},$

(16)

where $\Delta A^{\mathrm{ML}}$ is the free energy change determined by the particle insertion or element substitution method described in the previous subsection.

The correction of MLFF errors through this TI is generally considered necessary to verify the accuracy of the results. However, in the present application, all MLFFs were highly accurate, and this correction resulted in only a very small change, altering the calculated real potentials by an average of 50 meV.

In both the training runs for generating the MLFFs and the production runs for computing free energy changes, the liquid solvent was modeled using a unit cell with a side length of 12.4 Å, containing 64 water molecules. For the TI calculations, a single ion per unit cell was inserted, as illustrated in Fig. 2. The choice of the 64-water-molecule system was based on a previous study Jinnouchi_npjComputMater_2024 that showed this system size achieved sufficient convergence

All calculations were carried out using the Vienna Ab initio Simulation Package (VASP). Kresse_PRB_1996 ; Kresse_CMS_1996 For the MLFF models, I employed the kernel-based algorithm detailed in previous publications. Jinnouchi_PRB_2019 ; Jinnouchi_JPCL_2020 Similar to the pioneering methods, Behler_PRL_2007 ; Bartok_PRL_2010 the potential energy is expressed as a sum of atomic energies. Each atomic energy is represented as a weighted sum of kernel basis functions of the radial and angular descriptors. Jinnouchi_JCP_2020 The weights on these kernel basis functions are optimized to best reproduce the FP training datasets. These optimizations are performed within a Bayesian framework, which allows for accurate predictions of energies, forces, stress tensor components, and their uncertainties, thereby enabling efficient on-the-fly sampling of reference structures during FPMD simulations of the target systems.

To calculate the solvation free energy for each ion using the particle insertion method, a single MLFF was generated. The parameters for the descriptors and kernel basis function optimized in the previous study Jinnouchi_JPCL_2023 were utilized. Their values are listed in Supplementary Table S1. The MLFF was trained on liquid water on-the-fly during a 100-ps NVT-ensemble heating MD simulation, where the temperature was increased from 300 K to 500 K. Subsequently, the MLFF was trained on-the-fly during another 100-ps NVT-ensemble heating MD simulation of the solution containing the target ion. The MLFF underwent further training on-the-fly during TI calculations along the coupling constant $\lambda_{\mathrm{II}}$, which seamlessly transitions the soft repulsive Gaussian potential to the full interaction potential via Eq. (7). In this training run, a 6-point Gauss-Lobatto quadrature was employed. At each integration grid point, a 100-ps NVT-ensemble MD simulation was performed. A single MLFF for each ion was also generated using the element substitution method. The initial conditions before the on-the-fly training during the TI calculations were the same as those for the MLFFs for the particle insertion method. After these training runs, the MLFF was further trained during TI calculations along the coupling constant $\eta$. In this training run, a 5-point equispaced grid was used, and a 100-ps NVT-ensemble MD simulation was performed at each integration grid point. All MD simulations utilized the Nose-Hoover thermostat. Nose_JCP_1984 ; Hoover_PRA_1985 During the heating runs, the temperature was linearly ramped up from 300 K to 500 K. For efficient sampling, the mass of hydrogen atoms was increased to 4.0 amu, and the timestep of the MD simulations was increased to 2.0 fs.

The FP calculations were performed using plane wave basis sets with a cutoff energy of 520 eV and the projector augmented wave (PAW) method. Kresse_PRB_1999 The PAW atomic reference configurations were 1s¹ for H, 2s²2p⁴ for O, 2s²2p⁵ for F, 2p⁶3s¹ for Na, 3s²3p⁵ for Cl, 3s²3p⁶4s¹ for K, 4s²4p⁵ for Br, 4s²4p⁶5s¹ for Rb, 5s²5p⁵ for I, 5s²5p⁶6s¹ for Cs, 6s²6p⁵ for At, and 6s²6p⁶7s¹ for Fr. The $\Gamma$-point was used for Brillouin zone integration. In all calculations, the exchange-correlation interactions between electrons were modeled using the revised Perdew-Burke-Ernzerhof semi-local functional (RPBE+D3), Hammer_PRB_1999 augmented with Grimme’s dispersion interaction. Grimme_JCP_2010 ; Grimme_JCC_2011

In all on-the-fly sampling, a threshold spilling factor of 0.01 was used. If the spilling factor Miwa_PRB_2016 ; Jinnouchi_PRB_2019 exceeded 0.01, the MLFF prediction for the structure appearing in the MD simulation was considered to involve large errors. In such cases, an FP calculation was performed, the structure and the FP data were added to the training dataset, the MLFF was reconstructed, and the atomic positions were updated using the FP result. Otherwise, the MLFF prediction was used to update the atomic positions. This approach allowed us to bypass more than 99 % of FP calculations, enabling highly efficient sampling. Similar to the MLFFs in previous studies, Jinnouchi_PRB_2019 ; Zuo_JPCA_2020 ; Jinnouchi_JPCL_2020 ; Jinnouchi_JCP_2020 ; Jinnouchi_JCP_2021 ; Verdi_npjComputMater_2021 ; Lysogorskiy_npjComputMater ; Jinnouchi_JPCL_2023 ; Jinnouchi_npjComputMater_2024 ; Hijes_JCP_2024 the MLFFs achieve small RMSEs of 0.6 meV atom^-1 and 0.036 eV Å^-1 (see error distributions in Supplementary Figs. S1 to S4, enabling accurate and efficient computations of free energy changes.

The TI calculations by the particle insertion method were conducted using the variable transformation proposed by Dorner and co-workers: Dorner_PRL_2018

$\displaystyle\lambda_{\mathrm{I,II}}$

$\displaystyle=\left(\frac{x+1}{2}\right)^{\frac{1}{1-k}}.$

(17)

The variable $k$ was set to 0.5 in this study. After verifying the convergence with respect to the number of numerical integration grids [see examples in Fig. 3 (b)], the 12-point and 14-point Gauss-Lobatto quadratures were employed for the coupling constants $\lambda_{\mathrm{I}}$ and $\lambda_{\mathrm{II}}$, respectively. As in previous publications, the TI calculations were conducted in both forward and reverse directions along the coupling parameters to check the reversibility of the thermodynamic path. For the proton, which diffuses via the complex Grotthuss mechanism, 14 cycles of forward and reverse transitions were performed to achieve good statistics. For other ions, the number of cycles was set to 4-6.

The TI calculations using the element substitution method were conducted for the following pairs: H${}^{+}\leftrightarrow$Li⁺, Li${}^{+}\leftrightarrow$Na⁺, Na${}^{+}\leftrightarrow$K⁺, K${}^{+}\leftrightarrow$Rb⁺, Rb${}^{+}\leftrightarrow$Cs⁺, Cs${}^{+}\leftrightarrow$Fr⁺, F${}^{-}\leftrightarrow$Cl^-, Cl${}^{-}\leftrightarrow$Br^-, Br${}^{-}\leftrightarrow$I^-, and I${}^{-}\leftrightarrow$At^-. Similar to Eq.(17), the variable transformation from $\xi$ to $x$ was employed for the element substitution of H⁺$\leftrightarrow$Li⁺, where the integrand changed steeply near the endpoint of $\xi=1$. After checking convergence, 22-point grids were used for this transformation. For other element substitutions, where the integrands are very smooth, 6-point equispaced integration grids were used (see the convergence in Fig.3).

In all TI calculations, a 100-ps NVT-ensemble MD simulation was conducted at each integration grid point. The temperature was maintained at 300 K using the Nose-Hoover thermostat. The mass of hydrogen atoms was set to 2.0 amu, and the MD timestep was set to 1.0 fs. I ensured that the spilling factor for the structures appearing during the calculation remained below 0.02. This indicates that the structures encountered during the calculation are within the interpolation region of the training data.

The TI calculations to correct the results of the MLFFs via Eqs. (14) and (15) were performed using 3-point equispaced grids. At each integration grid point, a 10-ps NVT-ensemble MD simulation at 300 K was conducted. The same FP method used for the training runs was employed.

The local potential gaps $\Delta\phi$ for each solution system were computed from 1s levels of oxygen atoms in water molecules by using the method proposed in the previous study. Jinnouchi_npjComputMater_2024 ; Jinnouchi_arXiv_2024 The computed values of the 1s levels and $\Delta\phi$ are listed in Supplementary Table S2.

Figures 2 (a) and (b) show the real potentials of the proton (H⁺), alkali metal cations (Li⁺, Na⁺, K⁺, Rb⁺, Cs⁺, and Fr⁺), and halide anions (F^-, Cl^-, Br^-, I^-, and At^-) computed using the atom insertion method and the element transformation method. The values are also listed in Supplementary Table S3. Here, all the results shown in these figures are the FP results after applying the correction to the MLFF results using Eq. (16). As mentioned earlier, the difference between the MLFF and FP calculations was very small, with an average difference of 50 meV, which is also worth noting. Despite the completely different coupling pathways, the two methods yield identical real potentials within the statistical error bars. This result indicates good reliability and reproducibility of the computations.

The sum of the real potential for a cation and that for an anion equals the solvation free energy of a neutral ion pair, which can be precisely determined in experiments. Klots_JPC_1981 ; Tissandier_JPCA_1998 The calculated values are compared with the experimental ones in Table 1. While the calculations reproduce the experimental solvation free energies, close comparisons indicate that the RPBE+D3 functional systematically yields weaker ion-water interactions by an average of 0.1 eV per ion.

Although the results of the two TI calculations are identical, the element substitution method leads to a computationally more feasible thermodynamic pathway than particle insertion. As demonstrated by Fig.3 (a) (see also Ref. Jinnouchi_arXiv_2024 ), in the case of Na⁺ cation insertion, the TI calculation shows a steep change in the integrand near $\lambda_{\mathrm{I}}=0$, where the not-yet-interacting Na⁺ cation can come very close to other atoms. Additionally, the absolute value of the integrand is relatively large along the entire coupling pathway, and the integrand drops near $\lambda_{\mathrm{II}}=1$, where the fully interacting Na⁺ cation diffuses, escaping from the restraining model potential. In contrast, the element substitution from K⁺ to Na⁺ results in a very smooth and small integrand, as shown in Fig.3 (d). Consequently, the integrand converges more quickly with respect to the number of integration grids and the MD simulation time, as illustrated in Fig.3 (e) and (f). It is also noteworthy that the statistical error bars resulting from the element substitution are negligibly small. Therefore, once the real potential of one ion, such as K⁺ in this example, is computed using the particle insertion method, the real potential of another ion, Na⁺, can be easily obtained.

Now, I detail the solvation structures. Radial distribution functions (RDFs) ($g_{\mathrm{XO}}$) and running integration numbers (RINs) ($n_{\mathrm{XO}}$) between ions and oxygen atoms in water, computed by the MLFF and FP methods, are shown in Figs.4 and 5. Distances between the ion and the first- and second solvation shells, as well as the number of water molecules in these solvation shells, are listed in Supplementary Tables S4 and S5 for comparison with previous results. For all ions, the MLFFs reproduce the RDFs calculated by the FP method. The first peak in $g_{\mathrm{XO}}$ for both alkali metal cations and halogen anions shifts outward as the atomic number increases. The peak positions are in good agreement with previous results obtained by MD simulations using conventional force fields and FP methods, Dang_JCP_1993 ; Tuckerman_JCP_1995 ; Lee_JPC_1996 ; Toth_JCP_1996 ; Koneshan_JPCB_1998 ; Ramaniah_JCP_1999 ; White_JCP_2000 ; Raugei_JACS_2001 ; Lyubartsev_JCP_2001 ; Hribar_JACS_2002 ; Grossfield_JACS_2003 ; Ayala_JCP_2003 ; Izvekov_JCP_2005 ; Heuft_F_JCP_2005 ; Heuft_I_JCP_2005 ; Lamoureux_JPCB_2006 ; Dangelo_IC_2010 ; Bankura_JCP_2013 as well as experimental results. Lawrence_JCP_1967 ; Narten_JPC_1970 ; Caminiti_JCP_1980 ; Musinu_JCP_1984 ; Copestake_JP_1985 ; Skipper_JP_1989 ; Neilson_ARPC_1991 ; Cartailler_JP_1991 ; Ohtaki_CR_1993 ; Tanida_JCS_1994 ; DAngelo_JCP_1994 ; Enderby_CSR_1995 ; Fulton_JCP_1996 ; Beudert_JCP_1997 ; Wallen_JPCA_1997 ; Ferlat_PRB_2001 ; Kameda_BCSJ_2006 ; Filipponi_PRL_2003 ; Bowron_JP_2009 ; Mähler_IC_2012 However, close examination indicates that the RPBE+D3 functional appears to yield slightly larger solvation shells compared to those obtained by experiments and previous simulations using other functionals. This trend is consistent with the weaker ion-water interactions observed in the real potentials.

The Born model suggests that the solvation free energy is inversely proportional to the ionic radius:

$\displaystyle\Delta G_{\mathrm{solv}}$

$\displaystyle=-\frac{1}{2r}\left(1-\frac{1}{\epsilon_{\mathrm{r}}}\right),$

(18)

where $r$ is the ionic radius, and $\epsilon_{\mathrm{r}}$ is the relative permittivity of the solvent. This inverse dependence is observed in the real potentials of both cations and anions versus the position of the first solvation shell, as shown in Fig. 6. However, the simple model cannot explain why the real potential of an anion is lower than that of a cation of similar size, as shown in Fig. 6. In a previous study using a polarizable force field, Grossfield Grossfield_JCP_2004 discovered that the solvation of anions does not cause as much disruption in the surrounding water structure as the solvation of cations does. The same trend is observed in my MD simulations using the MLFFs. As shown in Fig. 7, the RDF between water molecules in the first solvation shell of a cation and other water molecules differs significantly from the RDF of bulk water without ion. Similar to the RDFs reported in previous studies, the height of the first peak diminishes to roughly the same level, and the first minimum disappears, especially for smaller ions. As pointed out by Grossfield, the reduction in the first peak is due to the volume excluded by the ion, and the diminished minimum indicates that water molecules in the first solvation shell have other water molecules packed around them without forming typical hydrogen bonds. The differences from bulk water gradually diminish with increasing ionic radius, indicating that the original water structure is restored as ion-water interactions decrease. In other words, these results suggest there is an energetic trade-off between stabilizing cation-water interactions and stabilizing water-water interactions. In contrast, the RDFs for halide anions are very similar to those of bulk water. Although the height of the first peak diminishes due to the volume excluded by the ion, the first minimum and second maximum appear at positions identical to those in bulk water. These results indicate that anion solvation does not significantly disrupt the water-hydrogen bonding network.

In summary, the MLFFs accurately replicate the dependence of the real potential of cations and anions on ion size and reproduce the differences in solvation structure between cations and anions, as discovered by a polarizable force field based on physical laws.