1 Introduction↩︎

Magnesium ions play a fundamental role in numerous biological processes, serving as stabilizing agents for proteins and RNA and as essential cofactors in a wide range of enzymatic reactions [1]–[6]. Classical molecular dynamics (MD) simulations have become a standard tool for investigating Mg\(^{2+}\) in complex aqueous and biological environments, and considerable effort has been devoted to the development and optimization of ion force fields to reproduce experimental data  [7]–[13]. Despite these efforts, classical force field simulations often fail to simultaneously reproduce key experimental structural and thermodynamic properties such as the hydration free energy and the structure of the first hydration shell [10], [11], [14], [15]. Additional limitations include unrealistically slow exchange kinetics in the first hydration shell of Mg\(^{2+}\) and excessive binding of divalent ions to nucleic acids [16], [17]. These limitations arise from the functional form of standard, fixed-charge classical force fields, which cannot explicitly capture quantum many-body effects such as polarization and charge transfer. Although approximate corrections, including modified combination rules [18], charge scaling approaches [19]–[21], and additional empirical terms in the pair potential [22], have been proposed, these methods do not provide a rigorous treatment of the underlying many-body effects.

Machine-learned interatomic potentials have emerged as a powerful alternative to conventional classical force fields [23]–[25]. In particular, neural network potentials (NNPs) provide a promising framework for extending near electronic-structure accuracy to the time and length scales required for the simulation of complex condensed-phase and biomolecular systems [24], [26]–[28]. A major breakthrough demonstrating the capabilities of NNPs has been achieved in simulations of liquid water, where NNP-based models trained on ab initio reference data were shown to reproduce structural, thermodynamic, and dynamical properties with near first-principles accuracy [23], [25], [29], [30]. In this way, NNPs can bridge the gap between electronic-structure calculations and experimentally relevant simulation scales, which remain inaccessible to conventional ab initio molecular dynamics (AIMD) due to its prohibitively high computational cost [31]–[33]. Recently, NNPs have also been developed for aqueous ionic systems and electrolyte solutions, including Mg\(^{2+}\), MgCl\(_2\), NaCl, Na\(_2\)SO\(_4\), and ZnCl\(_2\) solutions [34]–[38].

However, NNPs are not free of limitations, as they necessarily reflect the approximations and deficiencies of the underlying electronic-structure reference method used for training. For instance, the structural, thermodynamic, and dynamical properties of water and aqueous electrolytes are known to depend sensitively on the choice of exchange–correlation functional, the treatment of dispersion interactions, and subtle error cancellation effects. [36], [39]–[41]. Hence, careful validation against experimental reference data is essential. For classical force fields, parameterization and validation are typically based on experimental ion properties such as the solvation free energy, Mg\(^{2+}\)–oxygen distances in the first hydration shell, hydration numbers, activity derivatives, and water exchange kinetics in the first solvation shell. However, a comprehensive and systematic validation of NNPs for aqueous electrolyte solutions against a broad range of experimental observables has not yet been performed.

In this work, we present a robust MACE-based NNP for aqueous MgCl\(_{2}\) solutions and systematically investigate a broad range of structural, thermodynamic, and kinetic properties in comparison with established classical force fields and experimental reference data. To assess the sensitivity of the resulting properties to the underlying electronic-structure description, we train NNPs based on both the revPBE and hybrid revPBE0 exchange–correlation functionals [42]–[44] with D3 zero-damping dispersion corrections [45]–[47], motivated by the well-known functional dependence of water and electrolyte properties [36], [39], [40]. As a stringent test of whether NNPs can simultaneously reproduce the structure, thermodynamics, and dynamics of aqueous Mg\(^{2+}\) systems, we calculate solvation free energy, hydration structure, ion pairing, water-exchange rate and determine the mechanism of water exchange.

We find that both NNPs accurately reproduce the octahedral hydration shell structure of Mg\(^{2+}\) and predict a dissociative water-exchange mechanism, with computed exchange rate constants within one order of magnitude of experimental values. At the same time, the NNP-derived solvation free energies are significantly underestimated, indicating that despite the significant improvements offered by NNPs, important challenges remain in accurately describing long-range electrostatic interactions and charge-transfer effects in aqueous electrolyte systems.

2 Methods↩︎

2.1 Molecular Dynamics Simulations↩︎

Force Field Simulations: Initial classical force field molecular dynamics (MD) simulations were conducted using GROMACS 2023.05 [48] with the Mamatkulov-Schwierz \(\mathrm{Mg}^{2+}\)parameters [10] and the TIP3P flexible water model [49]. The simulation protocol comprised energy minimization, followed by 1 ns of NVT equilibration and 1 ns of NPT equilibration at 298.15 K and 1 bar. Final classical force field production runs were executed in the NPT ensemble for 5 ns. Detailed numerical parameters, including cutoff schemes, electrostatic treatments, and thermostat/barostat settings, are provided in the Supplementary Information.

Ab-Initio MD: All ab initio molecular dynamics simulations were performed using CP2K 2025 [50] within the Gaussian and Augmented Plane Wave (GAPW) framework [51]. Electronic structure calculations employed the revPBE exchange-correlation functional [42], [43] supplemented by Grimme’s DFT-D3 dispersion correction [45], [46] with zero damping, using TZV2P-GTH basis sets and GTH-PBE pseudopotentials [52], [53]. Self-consistent field convergence was achieved via the orbital transformation method with a DIIS minimizer.

The AIMD simulations were initiated from force field equilibrated structures following brief geometry optimizations. Molecular dynamics were conducted in the NVT ensemble at 298.15 K using the CSVR thermostat [54]. Detailed electronic and numerical parameters are provided in the Supplementary Information.

The initial NNP training dataset was produced by running five independent replicas for a pure water system (\(\mathrm{w}256\)) and a water system with solvated \(\mathrm{MgCl}_2\) (\(\mathrm{w}256\mathrm{MgCl}_2\)). For later testing of NNP transferability, simulations of a larger \(512\)-water system (\(\mathrm{w}512\mathrm{MgCl}_2\)) and each ion separately in the larger system (\(\mathrm{w}512\mathrm{Mg}\) and \(\mathrm{w}512\mathrm{Cl}\)) were performed. In order to obtain revPBE0 training data, we follow previous works by assuming that revPBE and revPBE0 phase spaces have significant overlap [55]. Under this assumption, we computed single-point revPBE0-D3/zd energies and forces for configurations extracted from revPBE-D3/zd AIMD trajectories at intervals of 5 fs, which is the stride we later use to subsample our training data.

NNP Simulations: All NNP simulations were conducted using OpenMM version \(8.2.0\) [56]. The neural network potential was implemented via the Message Passing Atomic Cluster Expansion (MACE) model [28]; parameter details are given in the next section. The Langevin equations of motion were integrated using the OVRVO integrator with time step rescaling [57]. The integration parameters were set to a time step of 1 fs and a friction constant of 1  ps. Simulations were typically performed in the NPT ensemble at a temperature of 298.15 K and a pressure of 1.0 bar, controlled by a Monte Carlo barostat. Enhanced sampling NNP-MD simulations, including metadynamics [58], [59] and umbrella sampling [60], were implemented via the PLUMED plugin [61], [62] to apply biasing potentials. Specific deviations from these standard simulation parameters, such as timestep adjustments or changes in ensemble conditions, are specified in the respective sections.

2.2 Neural Network Potential↩︎

Model Architecture and Training: We used the MACE architecture [28] to construct the NNP. MACE is an equivariant message passing neural network that uses higher-order body messages to describe atomic interactions. The model here consists of two message passing layers. The features are updated based on the local environment within a radial cutoff of 5.0 Å. The maximum symmetry order of the messages was set to \(L=0\), corresponding to an invariant model. We found no significant improvement setting \(L>0\) for our dataset. The number of feature channels was set to \(32\). These parameters were selected based on initial benchmarks showing a good balance between simulation performance and accuracy. The AIMD-derived \(\mathrm{w}256\) and \(\mathrm{w}256\mathrm{MgCl}_2\) datasets were split into training, validation, and test sets for both revPBE-D3/zd and revPBE0-D3/zd potentials, with isolated atomic energies (\(E_0\)) determined from DFT calculations. Models were trained separately on a single NVIDIA A100 GPU using the MACE reference implementation; detailed hyperparameters, batch sizes, and dataset splits are provided in the Supplementary Information.

Iterative Training: To generate additional training data, we performed NNP molecular dynamics simulations using the initially trained NNPs to explore phase space regions beyond the initial dataset, in particular transition states of the water exchange process. These included standard NPT simulations, high-temperature NVT runs, a 1 M \(\mathrm{MgCl}_2\) solution simulation, and two well-tempered metadynamics [58], [59] simulations biased along Mg\(^{2+}\)-oxygen and Mg\(^{2+}\)-Cl\(^-\) distances. We ran single-point DFT calculations on configurations extracted from these NNP-MD trajectories, added them to the initial dataset, and performed a second round of model training. Detailed simulation parameters are provided in the Supplementary Information.

2.3 Ion Properties↩︎

Structure of First Hydration Shell: To determine the first-shell Mg\(^{2+}\)-oxygen distance (\(R_1\)) and Mg\(^{2+}\) coordination number (\(n_1\)), we conducted ten independent equilibrium NNP-MD simulations for each system. These production runs were performed in the NPT ensemble at 298.15 K and 1.0 bar for a duration of 1 ns each. Other simulation parameters remained as stated above. The final configurations of these runs were then used as initial configurations for most of the following calculations.

Diffusion Coefficient: Diffusion coefficients (\(D_0\)) of the Mg\(^{2+}\) ion were computed from five additional NNP-MD equilibrium simulations performed in the NVT ensemble. These simulations utilized the same temperature (298.15 K) and time step (1 fs) as the NPT runs, lasting 1 ns per replica. Trajectories were analyzed using the Einstein relation to extract the self-diffusion coefficients. Specifically, we follow the workflow as explained in Grotz et al. [11].

Activity Derivative: Activity derivatives (\(a_{cc}\)) were determined from five NNP-MD equilibrium simulations of a 1.08 m (\(\approx \SI{1}{M}\)) \(\text{MgCl}_2\) solution in the NPT ensemble at 298.15 K and 1 bar. These production runs were performed for 1 ns using a simulation box with an edge length of 2.5 nm. While the integration timestep and temperature coupling were in line with the common NNP parameters stated above, the Monte Carlo barostat coupling was reduced to \(25\) steps to ensure proper density equilibration at the high ion concentration. Activity derivatives were determined as described in previous works [63].

Free Energy Calculations: Absolute solvation free energies (\(\Delta G_{\mathrm{solv}}\)) were computed using two distinct alchemical pathways to validate the NNP predictions. The first method followed an indirect thermodynamic cycle approach [64], where the free energy difference between the classical force field and the NNP level of theory was calculated. The second method employed a direct decoupling protocol based on neighborlist manipulation [65], allowing the solute interactions to be gradually switched off directly without transforming to a force field representation. The solvation free energy of neutral MgCl\(_2\) ion pairs \(\Delta G_{\mathrm{solv}}\) was obtained by first computing the single-ion solvation energies for Mg\(^{2+}\) and Cl\(^-\), applying corrections, and summing the results [11]. Detailed \(\lambda\)-window configurations, sampling durations, and correction schemes are provided in the Supplementary Information.

Exchange Mechanism from Transition Path Sampling: Flexible-length transition path sampling (TPS) simulations [66] were performed to characterize the mechanism of water exchange in the first hydration shell. The underlying molecular dynamics parameters followed the standard NNP-MD protocols outlined above. Initial reactive pathways were generated using steered molecular dynamics with a moving harmonic restraint on the Mg\(^{2+}\)-oxygen distance. The reaction progress was monitored by using an order parameter \(\xi\) to define two stable states: State A, where a water molecule resides in the first hydration shell, and State B, where an external water molecule has replaced it. Detailed parameters for the steering protocol, the exact mathematical definition of \(\xi\), and sampling settings are provided in the Supplementary Information.

Exchange Rate from Transition Interface Sampling: Transition interface sampling (TIS) simulations [67] were performed using the DFT-trained NNPs to calculate the rate constants for water exchange. Following the TPS section, we employed the same underlying NNP molecular dynamics parameters, order parameter \(\xi\), and stable state definitions. The flux through the first interface was estimated from unbiased MD trajectories. Initial reactive pathways for the interface ensembles were generated using steered molecular dynamics with harmonic restraints on two separate Mg\(^{2+}\)-oxygen distances and coordination numbers. Path sampling was conducted across multiple interfaces with parallel path swapping between ensembles to enhance sampling efficiency [68]. Detailed interface spacing, the functional form of the shooting weight function, and all numerical parameters are provided in the Supplementary Information.

Potential of Mean Force: The potential of mean force (PMF) along the Mg\(^{2+}\)–oxygen distance was computed using umbrella sampling [60] in PLUMED [61], [62] combined with the readily available NNP-MD equilibrium simulations. To reduce the computational load, umbrella windows were only run for regions missing in the PMF reconstructed from equilibrium data. The free energy profiles were reconstructed using the weighted histogram analysis method (WHAM) [69]. Detailed parameters regarding the window placement, spacing, and biasing force constants are provided in the Supporting Information.

Metadynamics: Well-tempered metadynamics [58], [59] was used to explore the free energy landscape along the Mg\(^{2+}\)–oxygen distance and the ion coordination number. The simulation parameters followed the standard NNP protocol outlined above and were performed using PLUMED [61], [62]. Detailed simulation and biasing parameters are provided in the Supporting Information.

Comparison to Classical Force Fields: For comparison of the NNP results, we employed two established non-polarizable Mg\(^{2+}\) force fields: the Mamatkulov-Schwierz Mg\(^{2+}\) parameters [10] and the microMg parameters [11] combined with the TIP3P water model. Both force fields were optimized to reproduce experimental solution properties including Mg\(^{2+}\)–oxygen distances, coordination number, activity derivative and solvation free energy. In addition these models provide a particularly useful comparison since microMg closely reproduces the experimental water exchange rate but follows an associative water exchange mechanism. The Mamatkulov-Schwierz model follows a dissociative exchange mechanism but significantly underestimates the exchange rate [16].

3 Results and Discussion↩︎

We first evaluate the accuracy and transferability of the DFT-trained NNPs. Subsequently, we present a systematic benchmark against experimental data and established classical force fields. The analysis covers Mg\(^{2+}\) hydration structure, self-diffusion, first-shell water exchange kinetics and mechanism, Mg\(^{2+}\)–Cl\(^{-}\) ion pairing, activity derivative, and solvation free energy (Fig. 1, Tab. ¿tbl:tab:Mg295properties?), thereby testing whether the NNPs can simultaneously describe structure, solvation and water-exchange properties in aqueous Mg\(^{2+}\) solutions.

3.1 Neural Network Potential Training↩︎

To obtain an accurate and transferable description of aqueous Mg\(^{2+}\), we develop a neural network potential (NNP) using the MACE framework. The training is performed in an iterative manner to progressively improve the robustness of the model and extend its coverage of relevant configurational space. An initial model is first trained only on configurations from unbiased AIMD simulations and is subsequently used to efficiently sample additional regions of phase space that are not represented in the original training data. The initial model showed good accuracy for equilibrium validation configurations. However, enhanced sampling using the preliminary NNP revealed significant errors in predicted energies and forces for configurations outside the configuration space regions represented in the initial training set (see Supporting Information Fig. S1 for a comparison before and after iterative refinement). This behavior is consistent with previous NNP studies, which showed that models can appear accurate while failing in rare-event or high-free-energy regions if such configurations are not represented in the training data [34], [35], [76], [77]. In the present case, the relevant rare event is water exchange in the first hydration shell of Mg\(^{2+}\), which occurs on the microsecond timescale and transiently gives rise to configurations with increased or decreased hydration numbers of the Mg\(^{2+}\) ion. Accurate representation of these configurations is essential for correctly describing both the exchange kinetics and the underlying exchange mechanism, as discussed in detail below.

After retraining on configurations generated through enhanced sampling with the preliminary NNP, the model showed a substantially improved description of transient water-exchange configurations with altered first-shell hydration numbers. Similar active-learning strategies were used by Juraskova et al. [34] and Ferretti et al. [35], where ligand-exchange or coordination-number sampling was explicitly included to improve the description of metal-solvent exchange regions. Importantly, the refined model also successfully extrapolates to larger and charged simulation boxes, particularly with respect to force predictions, as validated against DFT single-point calculations (see Supporting Information Fig. S1).

3.2 Structure of the First Hydration Shell and Self-diffusion Coefficient↩︎

The structure of the hydration shells is a key property of aqueous Mg\(^{2+}\). In particular, the strongly bound six water molecules in the first hydration shell govern many of the structural and kinetic properties of the ion. The distance to the oxygen atoms in the first hydration shell, \(R_1\), and the first-shell coordination number, \(n_1\), therefore provide direct measures of whether the model accurately reproduces the characteristic octahedral hydration structure.

Both the revPBE-D3/zd and revPBE0-D3/zd NNPs accurately reproduce the DFT (Fig. S2) and experimental first-shell properties of aqueous Mg\(^{2+}\) (Tab. ¿tbl:tab:Mg295properties?, Fig.  1). The predicted average first-shell distances, \(R_1\), of \(0.211\) nm for revPBE-D3/zd and \(0.206\) nm for revPBE0-D3/zd are both within the experimental uncertainty of \(0.209 \pm 0.004\) nm. In addition, both models yield a first-shell coordination number \(n_1 = 6\), consistent with previous simulations and experimental studies reporting an octahedral first hydration shell for the Mg\(^{2+}\) ion [31], [34], [35], [70], [78]. Overall, both DFT-trained NNPs accurately capture the key structural features of the Mg\(^{2+}\) first hydration shell.

The self-diffusion coefficient of Mg\(^{2+}\) provides a complementary measure of model accuracy, as it probes the mobility of the ion in solution. It therefore serves as an important test of whether a model that reproduces the correct hydration structure also yields realistic transport properties. The diffusion coefficient obtained with revPBE-D3/zd is in good agreement with experiment, yielding \(D_0 = 0.75 \times 10^{-5}\) cm\(^2\)/s compared to the experimental value of \(0.706 \times 10^{-5}\) cm\(^2\)/s [71]. In contrast, revPBE0-D3/zd predicts a lower diffusion coefficient of \(0.46 \times 10^{-5}\) cm\(^2\)/s, indicating significantly reduced Mg\(^{2+}\) mobility in solution. This trend is consistent with previous DFT studies of liquid water, which showed that revPBE-type GGAs tend to produce comparatively less structured and more mobile water, whereas revPBE0 does not necessarily improve the description of condensed-phase water dynamics despite the inclusion of exact exchange [36], [39].

Previous studies reported that revPBE-D3/zd reproduces bulk-water properties well, at least partly due to favorable error compensation, whereas revPBE0-D3/zd yields similar performance for bulk water and Cl\(^{-}\) hydration without systematically improving the description of cation–water interactions [36]. In addition, revPBE-D3/zd was shown to reproduce the experimental water density isobar more accurately than its hybrid counterpart [55]. In the context of ab initio Mg\(^{2+}\) simulations, Ferretti et al. [35] similarly noted that revPBE-D4 provides a reasonable description of aqueous Mg\(^{2+}\), although its performance depends sensitively on the treatment of dispersion interactions and associated error compensation. Overall, our diffusion results are consistent with the broader picture that revPBE-D3/zd yields favorable water-like dynamics, whereas revPBE0-D3/zd does not provide a systematic improvement for the structural and transport properties of Mg\(^{2+}\).

Table ¿tbl:tab:Mg295properties? and Fig. 1 also include a comparison to classical force fields, many of which reproduce the experimental Mg\(^{2+}\) diffusion coefficient with good accuracy. For such comparisons, however, it is important to note that diffusion coefficients obtained from classical force fields are commonly corrected for the viscosity of the underlying water model [79]. This procedure compensates for deviations in the water self-diffusion and is intended to isolate the contribution of the Mg\(^{2+}\) force field itself to the ion transport properties. In the present work, no viscosity correction was applied, such that the reported diffusion coefficients reflect the combined description of both Mg\(^{2+}\) and the underlying DFT-based water dynamics.

3.3 Water Exchange in the First Hydration Shell↩︎

Water exchange between the first and second hydration shells of metal ions is a fundamental kinetic process governing reactions in aqueous solution and biological systems. Reproducing both the rate and mechanism of water exchange represents a particularly stringent test for the accuracy of aqueous Mg\(^{2+}\) descriptions obtained with neural network potentials, since classical force fields have so far generally failed to reproduce both properties simultaneously [11], [16], [80].

In general, water exchange mechanisms are classified as associative or dissociative depending on whether the exchange proceeds through an intermediate or transition state with increased or decreased coordination number [80], [81]. Experimentally, the mechanism is inferred indirectly from the activation volume obtained from the pressure dependence of the exchange rate measured by \(^{17}\)O NMR spectroscopy, which indicates an interchange-dissociative (\(I_d\)) mechanism for Mg\(^{2+}\) ions  [73], [74].

A common first step for characterizing Mg\(^{2+}\) water exchange is the calculation of the potential of mean force (PMF) as a function of the Mg\(^{2+}\)–water oxygen distance \(r_\mathrm{Mg-O}\), which provides an initial estimate of the water-exchange free-energy barrier (Fig. 2).

The PMF obtained with revPBE-D3/zd exhibits a barrier height of \(12.6\,k_\mathrm{B}T\) (Fig. 2). In contrast, revPBE0-D3/zd predicts a substantially higher barrier of \(15.7\,k_\mathrm{B}T\). The revPBE-D3/zd barrier is in good agreement with results reported by Juraskova et al., who obtained a barrier of \(12.83\,k_\mathrm{B}T\) from umbrella sampling simulations  [34]. This agreement is notable because Juraskova et al. employed a cluster-trained MACE model based on a \(\omega\)B97X-D3BJ reference, whereas the present revPBE-D3/zd model was trained and applied for periodic MgCl\(_2\) solutions.

However, a one-dimensional PMF along \(r_\mathrm{Mg-O}\) alone is insufficient to resolve the underlying water-exchange mechanism [16]. To obtain a complete molecular picture of the water exchange, we therefore combine two-dimensional free energy landscapes along \(r_\mathrm{Mg-O}\) and \(n_1\) with transition path sampling and transition interface sampling simulations. Here, the choice of the oxygen atom in \(r_\mathrm{Mg-O}\) is arbitrary but remains fixed during the simulation.

The \(r_\mathrm{Mg-O}\)-\(n_1\) free energy landscapes (Fig. 3) indicate a dissociative pathway for revPBE-D3/zd and revPBE0-D3/zd due to the lower barrier at a reduced coordination number of \(n_1=5\). The TPS simulations confirm that the true reactive trajectories pass through a transient five-fold coordinated state (Fig. 3, top). This observation is in agreement with the interchange dissociative mechanism proposed by experiments [73], [75]. Based on the free energy landscape, the five-fold coordination state is lower in free energy for revPBE-D3/zd than for revPBE0-D3/zd. TPS shows a slightly longer-lived five-fold transition region for revPBE0-D3/zd which results in longer path lengths. These findings agree qualitatively with previous works using different functionals, which observed dissociation from the octahedral first shell before another water molecule entered during pulling [34], and that DFT-based potentials favor the five-fold state over heptacoordination [35].

We use TIS to obtain water-exchange rates based on true dynamical trajectories, thereby eliminating the need for any transition state theory (TST) approximation. For Mg\(^{2+}\) this is particularly important as TST along \(r_\mathrm{Mg-O}\) significantly overestimates the true exchange rate [16]. The TIS rate constants are \(82.8 \times 10^5\) s\(^{-1}\) for revPBE-D3/zd and \(1.29 \times 10^5\) s\(^{-1}\) for revPBE0-D3/zd, compared to the experimental values of (\(5.3 - 6.7) \times 10^5\) s\(^{-1}\) [73], [74]. Even though revPBE-D3/zd overestimates the experimental exchange rate and revPBE0-D3/zd slightly underestimates it, both results are within an order of magnitude, indicating good agreement.

In a previous study, Ferretti et al. [35] estimated rate constants from coordination-number free-energy barriers by exponential scaling relative to microMg, which was parametrized to reproduce the experimental exchange rate [35]. This approach results in approximately \(8.4 \times 10^7\) s\(^{-1}\) for revPBE-D4 and \(2.0 \times 10^6\) s\(^{-1}\) for revPBE-D\(^\mathrm{OPT}\). However, this approach based on TST assumes that the pre-exponential constant in the Arrhenius equation does not differ between force field and DFT descriptions, which may explain the deviation from our rates, since microMg water exchange is associative while the revPBE-D4/D\(^\mathrm{OPT}\) water exchange occurs through a dissociative mechanism.

In summary, the results show that water exchange in the first hydration shell of Mg\(^{2+}\) using NNPs with revPBE-D3/zd and revPBE0-D3/zd follows a dissociative mechanism, consistent with the experimentally assigned interchange-dissociative mechanism. Moreover, the calculated water-exchange rates are in reasonable agreement with experiments, representing an important step forward, as classical and polarizable force fields have generally failed to reproduce both properties simultaneously. [11], [16], [80]

3.4 Mg\(^{2+}\)–Cl\(^{-}\) Ion Pairing and Activity Derivatives↩︎

To assess ion pairing in MgCl\(_2\) solutions, we analyze the Mg\(^{2+}\)–Cl\(^{-}\) radial distribution function (RDF) and the activity derivative. The Mg\(^{2+}\)–Cl\(^{-}\) radial distribution function provides direct insight into local structure and the propensity for contact and solvent-shared ion pairs.

The Mg\(^{2+}\)–Cl\(^{-}\) RDF shows no stable inner-shell contact ion pair and no ion clustering under the studied conditions (Fig. 4). Instead, Mg\(^{2+}\)–Cl\(^{-}\) association occurs predominantly through solvent-mediated ion-pairs, as expected [82]. This result highlights the strong binding of first-shell waters to Mg\(^{2+}\) and indicates that Mg\(^{2+}\)–Cl\(^{-}\) association is primarily mediated by the hydration shell rather than direct contact.

Complementary to the Mg\(^{2+}\)–Cl\(^{-}\) RDF, the activity derivative provides a thermodynamic measure of ion association that can be directly compared with experiments. Moreover, it depends on the balance of ion–ion and ion–water interactions in solution. It hence provides a particularly stringent test of whether the NNP captures ion association with the correct balance between Mg\(^{2+}\)–Cl\(^{-}\) and hydration interactions.

The activity derivatives at salt concentration 1.08 m align well with experimental data (Tab. ¿tbl:tab:Mg295properties?, Fig.  1). Interestingly, the classical microMg and MS force fields with modified combination rules also reproduce the experimental activity derivative even though the propensity of Mg\(^{2+}\)–Cl\(^{-}\) solvent-shared ion pairs is significantly higher (Fig. 4). This effect is compensated by stronger hydration interactions, resulting in a similar overall thermodynamic response. The inclusion of many-body quantum effects therefore alters the balance between Mg\(^{2+}\)–Cl\(^{-}\) association and ion hydration relative to the classical force fields, while yielding a comparable activity derivative.

Still, it should be noted that the relative stability of contact and solvent-shared ion pairs can be highly sensitive to the electronic-structure reference, as shown for NaCl, and that revPBE-D3/zd and revPBE0-D3/zd can give similar ion-pair barriers while differing in structural details [36].

3.5 Solvation Free Energy↩︎

Solvation free energies provide a stringent thermodynamic test for ion models, since they depend on ion-water interactions, long-range electrostatics, polarization, and charge-transfer effects. They therefore constitute a common target in ion force field optimization. In contrast, solvation free energies have been explored much less for neural network descriptions, partly because absolute ion solvation free energies are challenging to compute for NNPs using alchemical transformation methods.

We calculate the solvation free energy \(\Delta G_{\mathrm{solv}}\) of neutral MgCl\(_2\) ion pairs using two different methods: one based on neighborlist decoupling and the other based on a thermodynamic cycle involving a transformation to a classical force field description (see Methods).

The NNP-derived \(\Delta G_{\mathrm{solv}}\) for revPBE-D3/zd and revPBE0-D3/zd shows significant deviations from experiment (Tab. ¿tbl:tab:Mg295properties?, Fig.  1), whereas the classical force fields reproduce the experimental values by construction, as solvation free energies were included among their parameterization targets [10]–[12], [78]. Specifically, the \(\Delta G_{\mathrm{solv}}\) is significantly smaller in magnitude, around \(-943\) to \(-984\) kJ/mol for revPBE-D3/zd and revPBE0-D3/zd, than the experimental value of about \(-2532\) kJ/mol [71]. The deviations from experiment may arise from several factors, including limitations of the underlying DFT reference, the absence of explicit long-range electrostatic interactions in the current local MACE description, and NNP artifacts associated with the net charge of the simulation box during alchemical transformations. In particular, the lack of explicit long-range electrostatics may limit the accuracy with which the long-range solvent response encoded in the DFT reference is reproduced and therefore contribute significantly to the observed deviations in ion hydration thermodynamics. This interpretation is consistent with recent discussions of long-range NNPs, which emphasize that local cutoff models can struggle with dilute ionic solutions, dielectric response, charge redistribution, and varying total charge states [83], [84].

The large deviation of the NNP-derived Mg\(^{2+}\) solvation free energies from experiment indicates a limitation of the present models for describing absolute ion hydration thermodynamics. At the same time, the NNPs accurately reproduce structural and dynamical properties, including hydration structure, diffusion, ion pairing, and water-exchange kinetics. These results suggest that ion hydration free energies constitute a particularly stringent benchmark for current NNPs. An important next step will be to assess whether recently developed long-range NNP architectures, such as MACE with charge equilibration and global charge states, can resolve this discrepancy while preserving the accurate local structure and dynamics obtained with the present models [36], [83]–[86].

4 Conclusion↩︎

Machine-learned interatomic potentials provide a promising route toward accurate simulations of ions in aqueous solution. For classical force field simulations, Mg\(^{2+}\) has remained a particularly challenging case, as force fields have failed in the past to simultaneously reproduce solvation free energy, size of the first hydration-shell, water-exchange rate, and exchange mechanisms, likely reflecting the absence of an explicit treatment of quantum many-body effects.

The aim of this work was therefore to train MACE NNPs on revPBE-D3/zd and revPBE0-D3/zd reference data and to systematically evaluate the ability of DFT-trained NNPs to reproduce key structural, dynamical and thermodynamic experimental properties of aqueous MgCl\(_2\) solutions. For this purpose, NNPs provide a unique opportunity to calculate experimentally accessible macroscopic solution properties which are inaccessible for computationally demanding DFT calculations. Moreover, the direct comparison with experimental data enables critical assessment of the underlying DFT reference and current limitations of NNP descriptions.

Both NNPs perform well with respect to the structure of the first hydration shell, the self-diffusion coefficient, and the activity derivative. The models reproduce the octahedral hydration shell structure and experimental Mg\(^{2+}\)-water oxygen distance, and revPBE-D3/zd yields a diffusion coefficient close to experiment. The slower diffusion obtained with revPBE0-D3/zd demonstrates that properties predicted by NNPs remain sensitive to the choice of reference density functional.

Transition path sampling and transition interface sampling provided unbiased mechanistic insight and approximation-free rate calculations for water exchange in the first hydration shell of Mg\(^{2+}\). This represents a remarkable achievement, as the combination of path-sampling techniques with NNPs extends DFT-level simulations to rare events occurring on the microsecond timescale, which remain inaccessible to conventional AIMD.

Both NNPs predict a dissociative water-exchange mechanism via a five-coordinated intermediate, in agreement with experimental findings [75]. Moreover, the calculated exchange rates are within one order of magnitude of experiment, representing a substantial improvement over classical force fields, for which the rate deviates by four orders in magnitude when reproducing the correct dissociative exchange mechanism [16].

Overall, the NNPs accurately reproduce a range of structural and dynamical properties of aqueous Mg\(^{2+}\). However, they significantly underestimate the experimental solvation free energy, revealing a clear limitation of the present local NNP for describing ion solvation thermodynamics, where long-range electrostatics are essential and net charged systems must be considered. The results hence support the use of DFT-trained NNPs for Mg\(^{2+}\) hydration structure, transport, and exchange kinetics, while demonstrating that accurate local structure and dynamics alone are insufficient to guarantee accurate hydration thermodynamics. Future work should therefore incorporate explicit long-range electrostatic treatment and charge-equilibration schemes and further assess the dependence of the results on the underlying electronic-structure reference.

Our results demonstrate that rigorous validation against experimental thermodynamic, structural, and kinetic solution properties is essential for assessing the predictive power of ion NNPs and establish a benchmark for the development of next-generation machine-learned potentials for aqueous electrolytes.

Data Availability↩︎

MACE models, simulation scripts, and TPS/TIS code are publicly available at https://git.rz.uni-augsburg.de/cbio-gitpub/Mg2-MACE.

5 Methods↩︎

5.1 Molecular Dynamics Simulations↩︎

5.1.1 Force Field Simulations↩︎

Classical molecular dynamics simulations for the production of initial ab-initio simulation configurations were performed using the GROMACS package SAbraham2015? with the Mamatkulov-Schwierz magnesium force field SMamatkulov2018? and the flexible TIP3P water model SJorgensen1983?. A flexible water model was chosen to better reproduce the behavior of water in ab initio simulations. Energy minimization was achieved via the steepest descent algorithm for a maximum of \(50000\) steps, with convergence criteria set at 100 kJ nm⁻¹. The simulations utilized the Verlet cutoff scheme with non-bonded interaction cutoffs of 1.0 nm, although if the simulation box was smaller, the cutoff was set to half the box length. Long-range electrostatic interactions were treated using the particle-mesh Ewald (PME) SDarden1993? method with a Fourier spacing of 0.12 nm.

A 1 ns NVT equilibration was performed at 298.15 K, followed by a 1 ns NPT equilibration at the same temperature and 1 bar pressure with a timestep of 1 fs. Temperature coupling during both stages utilized the stochastic velocity-rescaling thermostat SBussi2007? with a coupling time constant of \(\tau_\mathrm{T} = \SI{1.0}{\pico\second}\). For the NPT equilibration, the Berendsen barostat SBerendsen1984? was used for pressure coupling (\(\tau_\mathrm{P} = \SI{5.0}{\pico\second}\)).

Production runs were executed in the NPT ensemble for 5 ns. During these production trajectories, temperature was maintained at 298.15 K using velocity rescaling (\(\tau_\mathrm{T} = \SI{1.0}{\pico\second}\)), while pressure was controlled via the C-rescale barostat SBernetti2020? with a coupling time of \(\tau_\mathrm{P} = \SI{5.0}{\pico\second}\).

5.1.2 Ab-Initio MD↩︎

All ab-initio simulations were performed using CP2K version 2025 SHutter2014?. The GAPW method SLippert1999? was employed with a plane-wave cutoff of 700 Ry, five grid levels, and a relative cutoff of 50 Ry. Exchange-correlation effects were treated using the revPBE functional SPerdew1996?, SZhang1998? supplemented by Grimme’s DFT-D3 dispersion correction with long-range corrections and zero damping SGrimme2010?, SGrimme2011?. TZV2P-GTH basis sets combined with GTH-PBE pseudopotentials were used for all elements (q10 for \(\mathrm{Mg}^{2+}\), q1 for H, q7 for \(\mathrm{Cl}^{-}\), and q6 for O) SGoedecker1996?, SHartwigsen1998?. SCF convergence was achieved via the orbital transformation method with a DIIS minimizer, targeting an energy threshold of \(5.0 \times 10^{-7}\) Hartree.

Molecular dynamics were initiated from force-field equilibrated structures following a brief BFGS geometry optimization (five iterations). Simulations proceeded in the NVT ensemble at 298.15 K with a timestep of 1 fs. After an initial 100 steps minimization phase, equilibration utilized a strongly coupled CSVR thermostat SBussi2007? (\(\tau = \SI{5}{\femto\second}\)), followed by production runs with weaker coupling (\(\tau = \SI{100}{\femto\second}\)) lasting 1250 steps. System charge was adjusted according to the simulated ion species.

5.2 Neural Network Potential Training↩︎

For Neural Network Potential (NNP) training, the \(\mathrm{w}256\) and \(\mathrm{w}256\mathrm{MgCl}_2\) datasets were split into training, validation, and test sets for both revPBE and revPBE0 potentials. The training set contained \(2500\) configurations, the validation set contained \(125\) configurations (\(5\%\) of the training data), and the test set contained \(250\) configurations. Models were trained separately on a single NVIDIA A100 GPU for \(250\) epochs with a batch size of \(6\) and a learning rate of \(0.01\). Energy and force weights were set to \(\lambda_E = 1\) and \(\lambda_F = 1000\), and Stochastic Weight Averaging along with Exponential Moving Average of the weights were enabled.

5.2.1 Iterative Training↩︎

The following simulations were performed using the initial revPBE and revPBE0 neural network potentials to generate new configurations:

• NPT: 1 ns at 298.15 K and 1.0 bar for \(\mathrm{w}256\) and \(\mathrm{w}256\mathrm{MgCl}_2\) with stronger barostat coupling (every 25 steps).

• High temperature: 0.125 ns NVT at 500 K for \(\mathrm{w}256\) and \(\mathrm{w}256\mathrm{MgCl}_2\) with a 0.5 fs timestep.

• High concentration: 1 ns NVT for a 1 M \(\text{MgCl}_2\) solution.

• Metadynamics (\(r_\mathrm{Mg-O}\) + \(n_1\)): 2.5 ns NVT biased along the Mg–O distance and coordination number (bias factor 10, hill height 1.0 kJ mol⁻¹, stride 100 steps, adaptive width \(\tau = \SI{500}{\femto\second}\), \(\sigma_\mathrm{min}=\SI{0.02}{\nano\meter\squared}\) and 0.05).

• Metadynamics (Mg-Cl Distance): 2.5 ns NVT biased along the Mg–Cl distance and coordination number, using identical metadynamics parameters.

From each run, 50 evenly spaced configurations (100 for metadynamics) were extracted, their energies and forces computed via single-point DFT calculations, and added to the initial dataset. This augmented dataset was then used to perform a second round of model training using the same protocol as in the first round. Figure 5 shows the network’s accuracy in predicting energies and forces after the first training cycle, before the second cycle, and after the second cycle. The middle column illustrates clearly the initial lack of performance in unknown regions of configuration space explored via enhanced sampling, while the last column shows that accurate predictions can be recovered in a second training round. As a further validation, Figure 6 shows good agreement between the Potential of Mean Force (PMF) obtained via umbrella sampling with the NNP and the partial PMF calculated directly from distance histograms of the DFT simulation.

5.3 Ion Properties↩︎

5.3.1 Free Energy Calculations↩︎

Detailed numerical parameters for the solvation free energy calculations are provided here. For both alchemical approaches, equilibrium simulations were performed at discrete \(\lambda\)-windows using the NNP simulation parameters defined in the main text, such as temperature and pressure coupling. Each \(\lambda\)-window was simulated for 50 ps with a timestep of 1 fs. The following simulation protocols were repeated for w512Mg and w512Cl, based on which the solvation free energy of the Mg\(^{2+}\)-Cl\(^{-}\) ion pair was estimated.

In the thermodynamic cycle method SKarwounopoulos2024?, the correction term \(\Delta G_{\mathrm{MM} \to \mathrm{NNP}}\) was obtained via thermodynamic integration. Sampling was conducted across \(6\) equidistant windows with \(\lambda \in \{0.0, 0.2, 0.4, 0.6, 0.8, 1.0\}\), where \(\lambda=0\) represents the pure force field and \(\lambda=1\) the pure NNP potential. The force field side of the thermodynamic cycle (\(\Delta G_{\mathrm{MM}}\)) was set up in the same way as the initial force field simulation setup described in section S1.1 employing Mamatkulov-Schwierz Mg\(^{2+}\) with flexible TIP3P water. The term \(\Delta G_{\mathrm{MM}}\) was recalculated using the protocol described in Mamatkulov et al. SMamatkulov2018?, but did not show significant differences to the value reported for rigid TIP3P water. Correction terms as described in Grotz et al. SGrotz2021? were applied to \(\Delta G_{\mathrm{MM}}\). Therefore, no further corrections were applied to \(\Delta G_{\mathrm{MM} \to \mathrm{NNP}}\).

For the neighborlist decoupling method SPicha2025?, the solute-solvent interactions were scaled by modifying the interatomic distances provided to the NNP using a linear_to_cutoff shifting scheme. For this approach \(21\) \(\lambda\)-windows were utilized, distributed as \(\lambda \in \{0.0, 0.025, \dots, 0.25, 0.3, 0.35, 0.4, 0.5, \dots, 0.9, 1.0\}\), where \(\lambda=1\) represents the decoupled state. Free energy differences were subsequently estimated using the BAR method SBennett1976?. Corrections were applied as described in Grotz et al. SGrotz2021?.

5.3.2 Transition Path Sampling with NNP↩︎

Initial reactive pathways were generated using steered molecular dynamics by applying a moving harmonic restraint on the distance between an inner-shell oxygen atom and the magnesium ion. The bias center was shifted linearly from 0.23 nm to 0.38 nm over 100 ps with a force constant of 50000 kJ mol⁻¹ nm². We define an order parameter \(\xi = \exp(-10 r_2) - \exp(-10 r_1)\), where \(r_1\) denotes the distance between the magnesium ion and the pulled oxygen atom, and \(r_2\) represents the distance to the sixth nearest oxygen neighbor excluding the oxygen assigned to \(r_1\). Stable states A and B were defined by thresholds \(\xi < -0.08\) and \(\xi > 0.08\), respectively. Sampling was performed over 1250 trials with a maximum trajectory length of 10 ps, and shooting points were selected uniformly along the trajectories.

In addition to the associative/dissociative classification discussed in the main text, water exchange can also be categorized by the angle between the incoming and outgoing water molecules (Fig. 7). This results in two pathways, a direct mechanism (narrow angle) and an indirect mechanism (wide angle). Both pathways are observed across both DFT functionals, indicating that there is no clear preference for a specific pathway for the NNP-described magnesium ion.

5.3.3 Transition Interface Sampling with NNP↩︎

The flux through the first interface (\(\lambda_0 = -0.08\)) was determined from unbiased MD simulations with a duration of 0.05 ns per replica across \(5\) replicas. Initial reactive pathways were generated using steered molecular dynamics over 0.1 ns, applying harmonic restraints to two magnesium–oxygen distances \(r_1\) and \(r_2\), and a coordination number restraint (\(k = \SI{10000}{\kilo\joule\per\mol\nano\meter\squared}\), \(n_\mathrm{ref}=5.1\)). The bias centers for \(r_2\) and \(r_1\) were shifted from 0.21 nm to 0.5 nm and from 0.5 nm to 0.21 nm, respectively, over 100 ps with a force constant of 15000 kJ mol⁻¹ nm².

Path sampling was conducted across 11 interfaces ranging from \(\lambda_0 = -0.08\) to \(\lambda_{10} = 0.00\) with a spacing of \(0.008\). Each interface ensemble was sampled for \(1100\) trials with a maximum trajectory length of 10 ps. Shooting points were selected using a biased weight function SJung2017? \(w(\xi)\) defined by the interface position \(\lambda_i\), a capping value \(\xi_\mathrm{cap} = 0.02\), and a falloff parameter \(\alpha = 0.015\): \[\begin{align} w(\xi) = \begin{cases} 0 & \xi \geq \xi_\mathrm{cap} \\ 1 & \lambda_i \leq \xi < \xi_\mathrm{cap} \\ \exp\left(-\frac{|\xi - \lambda_i|}{\alpha}\right) & \xi < \lambda_i \end{cases} \end{align}\] Parallel path swapping was employed between adjacent interface ensembles after every trial to improve efficiency and reduce correlation SvanErp2007?.

5.3.4 Activity Derivatives with NNP↩︎

Activity derivatives were calculated from NPT radial distribution functions using Kirkwood–Buff theory SKirkwood1951?, as detailed for magnesium in a previous work by Grotz et al. SGrotz2021?. All ion-ion and ion-water radial distribution functions involved in this calculation are shown in figure 8. We follow the protocol described in Grotz et al. SGrotz2021?, only the normalization of the radial distribution functions prior to the numerical Kirkwood–Buff integral calculation is adjusted.

In general, Radial Distribution Functions (RDF) from molecular dynamics simulations contain finite size effects and noise due to finite simulation boxes and time. This can lead to the RDFs not converging to unity at the cutoff distance. Without correct normalization, the Kirkwood–Buff integrals are numerically unstable and very prone to diverge at larger distances. Therefore, the common approach is to scale the simulated RDFs \(g^\mathrm{sim}_{ij}\) using an estimated normalization factor \(f\) \[\begin{align} g_{ij}(r) = f g^\mathrm{sim}_{ij}(r)\, , \end{align}\] so that \(g_{ij}(r)\) is guaranteed to converge to unity at the end of the cutoff. The scaling factor \(f\) is usually obtained numerically from an average over the tail of the RDF.

For NNP simulations, it is particularly challenging to obtain low-noise RDFs due to the high computational cost of simulations. Instead of averaging to obtain \(f\) to improve convergence, we screened a grid of possible normalization values between \(0.95\) and \(1.05\) and picked the normalization value that minimizes the variance of the calculated Kirkwood–Buff integral after a cutoff of 1.62 nm. While this has close to no effect on force field-derived RDFs from long MD simulations, it leads to a much more stable normalization and hence convergence of the numerical integrals for noisy RDFs.

5.3.5 Potential of Mean Force with NNP↩︎

Umbrella sampling simulations were performed along the magnesium–oxygen distance to compute the PMF STorrie1977?. For the revPBE NNP, seven windows were employed with restraint centers positioned between 0.285 nm and 0.335 nm, spaced at intervals of 0.005 nm. For revPBE0, windows were spaced by 0.005 nm between 0.275 nm and 0.325 nm. Each window was biased using a harmonic restraint with a force constant of 50000 kJ mol⁻¹ nm². The umbrella sampling simulations were run for 2 ns per window in the NPT ensemble. The weighted histogram analysis method (WHAM) SKumar1992? was used to reconstruct the free energy profile from the combined equilibrium and umbrella sampling trajectory data. In most umbrella sampling analysis tools, this can be achieved by supplying the equilibrium data as a separate umbrella window with a zero force constant.

5.3.6 Metadynamics with NNP↩︎

Well-tempered metadynamics SLaio2002?, SBarducci2008? calculations were performed with a total duration of 7.5 ns. To constrain the sampling to the magnesium–oxygen distance of interest, a harmonic upper wall was applied to the \(r_\mathrm{Mg-O}\) coordinate at 0.7 nm with a force constant of 25000 kJ mol⁻¹ nm². The metadynamics bias was added by depositing Gaussian hills of initial height 0.8 kJ mol⁻¹ every 100 steps using a bias factor of \(10\). Fixed Gaussian widths of 0.015 nm were used for the distance \(r_\mathrm{Mg-O}\), and \(0.05\) for the coordination number. The first-shell coordination number \(n_1\) is evaluated using a rational switching function characterizing magnesium–oxygen contacts within a radial cutoff of 0.3 nm, utilizing switching exponents of \(n=12\) and \(m=24\), with all distance calculations accelerated via a neighbor list (\(N_{\mathrm{LIST}}\)) featuring a cutoff of 0.6 nm updated every 25 simulation steps.

References↩︎