Modelling ligand exchange in metal complexes with machine learning potentials

Metal ions are irreplaceable in many areas of chemistry, including (bio)catalysis, self-assembly and charge transfer processes. Yet, modelling their structural and dynamic properties in diverse chemical environments remains challenging for both force fields and ab initio methods. Here, we introduce a strategy to train machine learning potentials (MLPs) using MACE, an equivariant message-passing neural network, for metal-ligand complexes in explicit solvents. We explore the structure and ligand exchange dynamics of Mg 2 + in water and Pd 2 + in acetonitrile as two illustrative model systems. The trained potentials accurately reproduce equilibrium structures of the complexes in solution, including different coordination numbers and geometries. Furthermore, the MLPs can model structural changes between metal ions and ligands in the first coordination shell, and reproduce the free energy barriers for the corresponding ligand exchange. The strategy presented here provides a computationally efficient approach to model metal ions in solution, paving the way for modelling larger and more diverse metal complexes relevant to biomolecules and supramolecular assemblies.


Introduction
Metal ions have a central structural and functional role in many molecular systems, including catalysts, supramolecular assemblies, and biomolecules.Due to their relevance, much work has been done to investigate the structure, kinetics, and thermodynamic stability of metal complexes in solution, including the dynamics of metal-ligand exchange reactions. 1sing a variety of experimental techniques, including X-ray absorption spectroscopy, neutron scattering and nuclear magnetic resonance (NMR) spectroscopy, several mechanisms have been proposed to describe ligand exchange in the first coordination shell of the metal ion.These mechanisms range from dissociative (D), involving an intermediate of lower coordination number, to associative (A), proceeding through an intermediate of higher coordination number.However, these are extreme cases -in most instances, no such idealised intermediate exists, and instead, a concerted interchange mechanism with dissociative (I d ) or associate (I a ) characteristics occurs. 2,3f particular interest is ligand exchange with solvent, with metal aqua complexes being the most extensively studied. 4The a Chemistry Research Laboratory, University of Oxford, Oxford, OX1 3TA, UK.E-mail: fernanda.duartegonzalez@chem.ox.ac.uk b Center for Theoretical Chemistry, Ruhr University Bochum, D-44780 Bochum, Germany.E-mail: lars.schaefer@ruhr-uni-bochum.de ‡ Equal contribution to this work.† Electronic Supplementary Information (ESI) available: [details of any supplementary information available should be included here].See DOI: 00.0000/00000000.rate of this exchange depends on the nature of the metal ion, particularly ionic radii, charge, and coordination environment, ranging from 200 ps for Cs + to 300 years for Ir 3+ . 4Coordination with nonaqueous solvents such as alcohols, dimethyl sulfoxide (DMSO), acetonitrile (MeCN), and amides, has also been explored. 5mong the cations investigated, significant efforts have been made to study Mg 2+ complexes due to their prominent role in biology, including RNA folding, ATP hydrolysis, cellular signalling, and photosynthesis. 6][9] Water molecules in the first solvation shell are tightly bound to the cation and undergo exchange with the bulk solvent molecules on the microsecond timescale (k = 5.3 ×10 5 s −1 at 298 K) via a dissociative or interchange-dissociative mechanism. 10,116][17][18][19][20] Pd 2+ complexes have a square planar geometry defined by four coordinate bonds in equatorial positions.The complex can additionally interact with two more loosely bound ligands at the axial positions.In water, the [Pd(H 2 O) 4 ] 2+ complex has a Pd-O equatorial bond distance of 2.00-2.05Å, with a second solvation shell of 10 waters located between 4.02-4.40Å. 21,22 The Pd-O axial distance has been investigated by neutron diffraction 23 and extended Xray absorption fine structure experiments 24 and with different computational methods. 21,25,26The axial interaction distance is reported to range from 2.5 to 3.0 Å. Ligand exchange in Pd 2+ square planar complexes is suggested to occur via an associative mechanism involving a pentacoordinated trigonal bipyramidal transition state (TS), as suggested by Ligand Field Theory 27,28 and supported by DFT calculations. 295][36] Notably, the labile nature of Pd 2+ -ligand axial interactions is key for self-correction and optimal self-assembly. 37,38Pd 2+ -based metal-organic cages are commonly formed in MeCN solvent, 34,37,38 although water and DMSO are also widely-used.[Pd(MeCN) 4 ] 2+ has been characterised using single crystal X-ray diffraction, revealing a Pd-N bond length of 1.956 ± 0.008 Å. 39 NMR studies have explored MeCN ligand exchange, reporting reaction rates of k = 4.0 s −1 and k = 3.5 s −1 at 322 K. 40,41 Computational modelling of Mg 2+ and Pd 2+ cations has received significant attention, in particular for the former.7][48] MD simulations with non-polarizable FFs are the most widely applied option since they balance the cost and accuracy of the resulting dynamics.Here, the metal ions are modelled as single or a small set of point charges with the electrostatic, dispersion, and excluded volume interactions taken into account by a pairwise interaction potential.0][51] A 12-6-4 LJ potential has been developed to partially account for charge-induced dipole interactions via the r −4 term. 52,53However, none of the available Mg 2+ FF parameters can simultaneously reproduce all properties with sufficient accuracy.Moreover, given the focus on metal-water properties alone, these models cannot describe orbital-specific and anisotropic features important in many metal-containing protein or synthetic catalyst active sites, or even the properties of simple electrolytes. 54Polarizable FFs are in principle able to remediate the limitations of non-polarizable FFs, but they are less frequently used due to their time-consuming parameterization and increased computational cost, especially for exploring long-timescale processes. 42,50acKerell et al. 55 modified Mg 2+ parameters to describe its interaction with water, Cl -ions, and nucleic acids using a polarizable FF based on the Drude oscillator model.This approach uses QMcomputed interaction energies and geometries of hydrated complexes as reference as well as condensed-phase osmotic pressure calculations.Mg 2+ parameters for the AMOEBA force field were reported by Jiao et al. 56 and further refined by Piquemal et al. 57 However, both implementations experienced rapid water dissociation.This issue was addressed by Kurnikov and Kurnikova, 58 who introduced a distance-dependent polarization response for water.
The effect of various FFs on the ligand exchange mechanism in Mg[(H 2 O) 6 ] 2+ complex was extensively studied by Schwierz et al., using transition path sampling. 59,60They demonstrated that while the commonly used non-polarizable FFs correctly predict the dissociative characteristics of the mechanism of the water exchange (I d ), they tend to overestimate the free energy barrier, leading to a significantly slower reaction rate. 59In comparison, the polarizable FF Amoeba and specialized non-polarizable FF microMg lead to a preference for an associative mechanism (I a ) with a reaction barrier too low, leading to significantly faster reaction rates, further illustrating the complexity of the ligand exchange process. 60he interactions of Pd 2+ with water molecules have also been studied computationally with classical FF approaches.Sanchez Marcos et al. 21investigated the [Pd(H 2 O) 4 ] 2+ complex in water using MD simulations.They developed two intermolecular potentials to describe the interactions between Pd 2+ and the water molecules, one for the first solvation shell, fitted to interaction energies computed at the MP2 level on the gas-phase complex, and another for the hydrated ion-bulk water interactions by incorporating a continuum polarizable model to account for solvation effects. 21They suggested the presence of solvent molecules in the axial position located between 2.5 and 3.0 Å, referred to as a 'mesoshell'.The concept of the mesoshell has sparked debate within the scientific community, with recent studies suggesting that the structure of Pd 2+ aqua complexes in water should be interpreted under the 'extended first shell' paradigm. 25,44For example, utilising QM/MM methods, Adnan Ali Shah et al. 25,44 identified a weakly bound axial ligand (Pd-O distance of 2.8 Å), resulting in a broad peak in the RDF between the first and second solvation shells.Contrasting findings were reported by Chen et al. 26 using subsystem DFT AIMD simulations of the same complex.Their results indicated that water molecules rarely occupied the axial region.Instead, solvent molecules formed a protective "dome" on both sides of the square planar complex via strong hydrogen bonds, preventing the penetration of single water molecules from the axial directions.These studies provide alternative interpretations of the experimental data, underscoring the complex nature of axial interactions in Pd 2+ aqua complexes.
2][63] However, its high computational cost limits its use to small systems and picosecond timescale processes, often insufficient to obtain converged free energies and model ligandexchange processes.Machine learning potentials (MLPs) have emerged as promising alternatives to AIMD, reproducing accurate energies and forces from electronic structure reference calculations at a much lower cost. 646 and in electrolytes confined to nanoscale pores. 77aditionally, MLPs are trained using AIMD reference data under periodic boundary conditions (PBC).This approach inherently captures long-range interactions but incurs a high computational cost due to the large size of the system, primarily consisting of solvents.Consequently, this also limits the use of high levels of theory and restricts the achievable sampling.Previous works by our group 78,79 and others 80,81 have demonstrated the efficiency and accuracy achievable using cluster data for training.When combined with active learning (AL), which iteratively builds the training set based on a preliminary version of the trained MLP, this approach yields accurate and data-efficient MLPs at a low computational cost.In this study, we expand this protocol to model metal complexes in solution, using clusters of solvated metal ions for training.Specifically, we apply Atomic Cluster Expansion (ACE) 82 and its message-passing neu-ral network-based variant (MACE) 83 to two model systems, Mg 2+ complexes in aqueous solution, representing a strongly interacting and biologically relevant metal ion, and Pd 2+ complexes in acetonitrile (MeCN), a transition metal relevant for supramolecular chemistry in non-aqueous solvents. 34,37,38Ligand exchange in these complexes proceeds via different mechanisms, allowing us to investigate the capability of the MACE potentials to model structural and energetic features characteristic of both processes.

ACE and MACE machine-learning potentials
In this work, MLPs were trained using linear regression with the ACE 84 descriptor and its variant MACE, in which ACE is combined with an equivariant message-passing neural network architecture. 83he ACE descriptor 85 builds on a traditional many-body expansion, where the Potential Energy Surface (PES) of the system is expressed as a sum of different body-order interactions, including two-body, three-body, and higher-order interactions depending on the truncation.Although this approach is physically motivated, it is limited to modestly-sized molecular systems, as the computational cost of evaluating the energy scales exponentially with system size, making it impractical to consider inter-actions beyond the three-body order.ACE overcomes this limitation by introducing the concept of atomic neighbour density, where the energy of each atom depends on the many-body interactions with its N neighbours within a defined cut-off radius.The validity of this concept is based on the assumption that the energy of each atom only depends on its local environment. 75,86econd, it projects these densities onto physically invariant basis functions.This procedure ensures that the evaluation cost of many-body terms scales linearly, rather than exponentially, with the number of neighbours, regardless of the body order.A detailed description of the method is provided in Ref. 85.The capability of the ACE descriptor to accurately map the PES enables the use of simple linear regression for fitting, resulting in an accurate and data-efficient approach to train MLPs. 84ACE combines the ACE descriptor with an equivariant message-passing neural network architecture, 83 incorporating body-order contributions as node features.Using graph neural networks (GNNs), the body-order term and cluster region implicitly expand with the number of message-passing layers, resulting in a more accurate representation of atomic environments. 838][89] Linear ACE provides high accuracy in low-data regimes, making it particularly suitable for use in the early stages of AL, where typically small data sets are used. 78,79In this work, we use linear ACE to build the training data sets using the AL loop, while MACE is used to expand the data sets and fine-tune the final potentials used for production MD simulations.

Workflow
The workflow presented here builds on our previous work on automated AL strategy for modelling chemical reactions in explicit solvents (Fig. 1). 79The AL cycle is initiated from a small training set of approximately 10 structures.These data consist of gasphase molecules generated by random displacement from a QMoptimized structure or solvated clusters, obtained from MD simulations or random placement of molecules in a box.The structures are labelled with energies and forces computed at the reference level of theory.The initial training data set is extended using AL as follows: A first version of the MLP is generated and used to propagate several independent MD simulations, typically ten, for n 3 + 2 fs, where n is the index of the MD run in the AL loop, starting from 0. From these trajectories, new structures are selected using a similarity selector, 79 which identifies new structures to be included in the training based on the similarity between a global Smooth Overlap of Atomic Positions (SOAP) representation of data point p and configurations p ′ in the existing training set. 90The similarity vector, K, is defined as follows: where the components of the vector k(p • p ′ ) correspond to the SOAP kernel functions computed between the SOAP representation of new structure p 0 and the i-th configuration in the existing training data p i .Parameter ζ is a positive integer that increases the sensitivity of the kernel to changes in atomic position. 90The selector adds structures to the training set if the maximum value of their similarity vector, K, is smaller than the given threshold k T , i.e., max(K) < k T .The new structures are then labelled by the reference energy and forces, added to the training set and potential is retrained.If no structures are selected from the trajectories, the index n increases by one and a longer MD simulation is performed with the same potential.The AL procedure is repeated until it either reaches the maximum number of AL cycles or when no new structure is selected within the maximum AL time.Details on electronic structure and MD protocols are provided in the Computational Details section.

Data set preparation 2.3.1 Mg 2+ in water
The training data set for Mg 2+ consists of 344 structures.To increase the accuracy of the resulting potentials across systems with different sizes, the final dataset combines several subsets with different compositions, corresponding to three subsets: (i) [Mg(H 2 O) 6 ] 2+ complex in the gas phase (30 structures), (ii) Mg 2+ solvated in 45 to 53 water molecules in a spherical cluster with radius 7 Å (165 structures), and (iii) water clusters containing 45 to 49 water molecules placed in a spherical cluster with radius 7 Å (149 structures).The schematic representations of the structures are depicted in Fig. 1a.The cluster size was selected to be larger than the distance cut-off of the descriptors used in MLP, i.e., 6 Å, needed to cover the Mg-O distance range sampled in the dissociative mechanism.All datasets were trained using the energies and forces computed at ωB97X-D3BJ/def2-TZVP 91,92 level of theory as ground truth, which provides accurate estimation for structural and thermodynamic properties of large systems.The ACE MLP was used to generate the structures during AL unless specified otherwise.
[Mg(H 2 O) 6 ] 2+ complex in the gas phase AL was initiated from 10 structures obtained by a random displacement of [Mg(H 2 O) 6 ] 2+ complex in the gas-phase.The new structures were selected using the similarity selector with a SOAP cut-off of 5 Å and threshold of 0.999, with the maximum time in the active training loop set to 3 ps.This procedure led to the selection of 30 structures.

Water cluster subset
The initial structures of the bulk water system were prepared by classical MD simulation using TIP4P-Ew FF. 93 A cubic box (L= 15 Å) was solvated with 112 water molecules, minimized and equilibrated in an NPT ensemble (300 K and 1.0 bar) for 1 ns using Langevin dynamics and Berendsen barostat as implemented in the sander module of Amber-tools23. 94For the initial training set, 10 clusters of 7 Å radius containing 45 to 49 water molecules were cut from the equilibrated trajectory and labelled with the reference energies and forces (vide infra).The training set was then enhanced by AL using the similarity selector with a threshold of 0.9998 to avoid selection of too distorted structures, and a maximum time set to 10 ps to accommodate for water relaxation.This approach yielded an overall 149 structures.

Faraday Discussions Accepted Manuscript
Open Mg 2+ complex in a water cluster The initial structures for the training were generated using classical MD simulations including a Mg 2+ ion solvated by 325 water molecules in a box of 21.5 Å.The box was minimized and equilibrated to experimental water density using the same procedure as pure water.The system was modelled using the TIP4P-Ew water model combined with Li/Merz ion parameters in TIP4P-Ew water (12-6 HFE set). 5210 structures were randomly extracted from the equilibration trajectory and cut into a 7 Å sphere with a centre in the Mg 2+ cation.
The three subsets were combined to form the full data set used to initiate AL, which was performed in four phases.In the first phase, the ACE potential was trained on the full data set and used in AL with MLP-MD initiated from Mg 2+ -water cluster (structure selected from a subset (ii)).During the MLP-MD dynamics step, the cluster was constrained in its spherical shape by a flat-bottom spherical bias potential set at a distance of 8 Å from the Mg 2+ cation, applying a harmonic restraint of 500 kcal/mol Å −2 .MD simulations during AL were performed in an NVT ensemble at 400 K. New structures were selected using the similarity selector with a threshold of 0.999, with maximum time in AL set to 7 ps.This longer AL time compared to the gas phase complex was used to account for the higher flexibility of water molecules outside the first solvation shell.The first phase collected 11 structures, which were added to the data set.
In the second phase, the AL was initiated using one of the clusters generated during the previous phase, with the radius reduced by placing the harmonic bias at a distance of 7.5 Å to increase the density of the solvent around Mg 2+ and enhance the sampling of the repulsive region of the potential.To avoid selecting overly distorted structures resulting from potentially unstable dynamics, the selector threshold was tightened to 0.9999.Simultaneously, the AL time was increased to 20 ps to further capture longer water dynamics.The second phase resulted in generating an additional 55 structures.
After these two phases of training, we tested the accuracy and stability of the resulting potential by comparing the MLP energies and forces with the ground truth data (for details see SI §S2) and conservation of the total energy in NVE dynamics.Despite the high accuracy of the trained ACE potential, as evidenced by a mean absolute deviation (MAD) of 0.79 meV/atom for energy and 46 meV Å −1 for forces (Fig. S1), the NVE and NVT MD simulations using PBC with the ACE potential were unstable.This included the formation of bubbles, followed by the system's collapse.These instabilities could be alleviated by introducing more radial functions into the descriptor and tuning the hyperparameters.However, we decided to change the model in the following phases from ACE to MACE, which is computationally more efficient.Indeed, the MACE potential provided stable NVE dynamics under PBC using the same training data set.Interestingly, the long NVE dynamics with the MACE potential promoted proton transfer of a water molecule in the first solvation shell, leading to the formation of [Mg(H 2 O) 5 (OH)] + and H 3 O + species.As these structures were underrepresented in the previous AL loops, resulting in larger prediction errors (see Fig. S2), we manually selected 35 structures along the NVE trajectory, cut them into 7 Å clusters, which contained the species, and added them to the training set for the third training phase.Apart from adding these data, we further repeated the AL loop to account for possible differences between the conformational space sampled by ACE and MACE potentials.The re-trained MACE potential was again found to provide an accurate estimate of energies and forces, 0.69 meV/atom and 29 meVÅ −1 (see Fig. S3).The preliminary NVT dynamics with this potential, however, led to fast dissociation of one H 2 O molecules from the first solvation shell, without no exchange with the bulk solvent.To correct for this behaviour and to ensure accurate exchange of the water molecules around the Mg 2+ , we completed the training set by adding 23 structures with H 2 O molecule dissociated to a distance above 3.0 Å (see Fig. S4) from Mg 2+ within the fourth training phase.

Pd 2+ in MeCN
MLP for Pd 2+ complex was trained using a total of 581 data points from the following subsets (Fig. 1b).As in the Mg 2+ case, the size of all clusters was selected in a way that the resulting cluster radius exceeds the 6 Å cut-off used in ACE and MACE descriptors.Unless specified otherwise, ACE was used as the ML model in all training phases, employing the ground truth potential TPSS0-D3BJ/def2-TZVP, 92,95,96 since this functional has shown good performance in reactions involving late-transition metals. 95MD simulations in the AL loops were performed in an NVT ensemble at 300 K. MD simulations longer than 1 ps used a flat-bottom spherical harmonic bias potential set at a distance from the Pd 2+ cation, applying a harmonic restraint of 100 kcal/mol Å −2 .The value of the harmonic restraint has been chosen to maintain the integrity of the cluster without creating artefacts from the pulling force.The onset distance of the flat-bottom potential was varied according to the size of the clusters so that the density of the cluster was close to the experimental density of the (bulk) liquid.
[Pd(MeCN) 4 ] 2+ •MeCN complex in the gas phase The transition state (TS) structure corresponding to pentacoordinated trigonal bipyramid was obtained at the TPSS0-D3BJ/def2-TZVP level of theory (Fig. 1b) and was used as a starting point for a 2D relaxed PES scan along the two Pd-N bonds involved in the ligand exchange.These bonds are of equal length in the TS (2.28 Å); while the length of the remaining three Pd-N bonds is equal to 1.95 Å.The PES scan resulted in a total of 25 structures (5x5 grid).In the following paragraphs, TS refers to the transition state structure of [Pd(MeCN) 4 ] 2+ •MeCN obtained from the PES scan mentioned above, while reactant state (RS) refers to the structure obtained from the geometry optimization of TS in gasphase with TPSS0-D3BJ/def2-TZVP.
[Pd(MeCN) 4 ] 2+ complex solvated in clusters of 60 MeCN molecules The initial structure for the subset was generated with the Quantum Cluster Growth (QCG) method 97 (as imple- mented in xTB (version 6.5.0) 98 and the lowest energy conformer was selected with CREST. 99.Then, the conformational space was explored via metadynamics simulations with xTB. 98,100Three simulations of 25 ps each were run in the NVT ensemble at 600 K, using the Cartesian root-mean-square-deviation (RMSD), with respect to a list of reference structures updated every 2 ps, as a collective variable. 100During the simulations, a flat-bottom spherical bias potential was applied to keep the cluster at the experimental density of liquid acetonitrile at room temperature. 100For each of the three trajectories, the first 600 fs were excluded and the remaining frames were merged in a single trajectory.From this, frames have been extracted every 240 fs resulting in 305 structures.

[Pd(MeCN) 4 ] 2+ •MeCN complex in clusters of 20 -30 MeCN molecules
The structures used as starting conformations in the multi-step AL approach described in this section have been generated with CREST.The TS and RS structures obtained in the gas phase (from the dataset (i)) were solvated by 20 to 30 MeCN molecules using the QCG method 97 and the lowest energy conformation of the cluster was selected with CREST, 99 as mentioned above.Different conformations were generated for the different AL training phases as follows.Firstly, the data set (i) was used as an initial data set in the AL loop, starting from TS solvated by 20 MeCN molecules; new structures were selected using the similarity selector with a threshold of 0.9999, with a maximum MD simulation time set to 750 fs, the time found to be required for transitioning from trigonal bipyramidal to square-pyramidal coordination geometries, observed from previous trials.This AL phase yielded 56 structures.
In the second phase, dataset (i) and the 56 structures obtained from the first phase were employed as the starting training set.The TS solvated by 30 MeCN molecules was used as a starting conformation for the new AL loop, with an MD simulation time of 1.5 ps and a SOAP selector threshold of 0.999.The choices to extend the MD simulation time and the number of solvating MeCN molecules were made to capture how the relaxation from the transition state would evolve on longer timescales and in the context of a larger solvation environment.The selector threshold was relaxed to avoid the selection of conformations too similar to the starting one since the fluctuations along the TS relaxation have been explored in the previous training phase.In this phase, 36 new structures were generated, expanding the dataset to 117 structures.
In the following phase, the RS solvated by 20 MeCN molecules was used as a starting structure, with the AL MD time extended to 5 ps to allow sampling of more distant regions of the PES.Therefore, the number of solvent molecules was decreased to lower the computational cost.30 new structures were generated and added to the previous data, yielding a total of 147 data points.Preliminary validation in the NVT ensemble under PBC showed artefacts in the description of the average structure of the system, in particular, the formation of void regions in the axial positions of the complex (Fig. S5).In line with the Mg 2+ case, we, therefore, decided to adopt the MACE for AL due to its computational efficiency and accuracy. 89urthermore, to enhance the stability of the potential during AL, we selected 150 structures from the dataset (ii) and merged them with the previous 147 structures, obtaining an extended starting dataset.A final AL loop was started from RS solvated by 20 MeCN molecules, using a similarity selector threshold of 0.9999 and simulation time set to 20 ps to ensure the stability of the potential on longer timescales.This training phase resulted in 129 new structures.Eventually, not considering the 25 structures from dataset (i), the multi-step AL approach described above yielded a total of 251 structures (147 + 129 -25).

Validation of the MLP
As a first step, we validated the accuracy and stability of the generated MACE potentials for both Mg 2+ in aqueous solution and Pd 2+ in MeCN.To evaluate the prediction error on unseen data, we generated an ensemble of testing structures by an MLP-MD simulation of the metal ion in a spherical cluster of solvent molecules.We then selected frames along the trajectory and performed a point-to-point comparison between the energies and forces computed at the ground truth DFT level of theory and MACE.
The testing data set for Mg 2+ cation in aqueous solution consisted of 51 structures collected over 50 ps NVT dynamics of Mg 2+ solvated in a cluster of 51 water molecules.The spherical shape of the cluster was kept by a harmonic spherical potential placed at 7.5 Å from Mg 2+ .Validation results are depicted in Fig. 2a.MACE demonstrates excellent performance in energies and forces, with MAD of 0.31 meV/atom and 18 meV Å −1 for energies and forces, respectively.
For Pd 2+ in acetonitrile, the MACE potential was tested on structures generated from 100 ps NVT dynamics using a cluster containing Pd 2+ and 30 MeCN molecules.The solvent molecules in the cluster were confined by a flat-bottom spherical bias potential with a radius of 10.0 Å, centred on Pd 2+ .From the trajectory, 51 structures were extracted.The MAD is 1.04 meV/atom and 8 meV Å −1 for energy and forces respectively.The good correlation with respect to the ground truth energies and the accuracy in forces suggest that the MACE potential closely reproduces the shape of the reference PES, with the higher energy MAD likely arising from a systematic shift in the absolute values, which has been previously reported in some instances with MACE. 89Overall, the MACE potentials provide a reliable prediction of energies and forces for both tested systems.
To further assess the performance of the MACE potentials in larger systems, we performed 100 ps MD simulations under PBC in the NVE ensemble for a system consisting of Mg 2+ with 145 water molecules in a 16.3 Å box and Pd 2+ with 159 MeCN molecules in a 24.0 Å box.In both cases, the MACE potential conserved energy, confirming the stability of the dynamics under PBC and on simulation times longer than the active learning time (Figs.S7  and S8).

Structural properties of the metal solvation shells
We evaluated the structural properties of the metal environment by computing the radial distribution functions (RDFs) between the metal ions and the coordinating solvent atom.For [Mg(H 2 O) 6 ] 2+ , the computed Mg-O RDF from 500 ps MD simulations (Fig. 3a) shows a first peak at 2.08 Å, in agreement with the 2.10 Å reported from X-ray diffraction and neutron scattering experiments. 7,101Integration of this curve results in a coordination number of 6, confirming the octahedral arrangement of this complex.A second, less well-defined peak is evident around r = 4.15 Å, corresponding to the second solvation shell.Integration of this peak results in a coordination number of 12, in agreement with experiments 7 and reported ab initio computations. 9In line with the known lifetime of the first solvation shell, the octahedral complex remained stable during the simulation time and no ligand exchange with the bulk solvent was observed.
As expected from the ligand field theory, 27 Pd 2+ forms a square planar complex with 4 MeCN molecules, with a Pd−N distance of 1.97 Å (Fig. 3b).This is in excellent agreement with the value obtained from the single crystal X-ray diffraction (SC-XRD) of the complex (Pd-N bond 1.96±0.01Å). 39 Two peaks around 3.3 Å, with a shoulder starting from 2.5 Å, and around 4.6 Å, are also observed.The latter peak is associated with the second solvation shell, with a coordination number of 8.The former peak corresponds to the interactions of acetonitrile molecules in the axial position, with a coordination number increasing from 4 to 6.The previous studies on Pd 2+ aqua complexes provided a foundation for understanding the axial interactions of Pd 2+ with MeCN.As mentioned in the introduction, two paradigms exist in the literature.In the mesoshell paradigm, the two axial ligands are symmetrically bound, resulting in a sharp peak between the first and second solvation shells. 21,23,24Conversely, the "extended first solvation shell" concept suggests more weakly bound axial ligands, leading to a presence of broad peak. 25,44The structural features observed in the RDF in Fig. 3b indicate that for Pd 2+ in acetonitrile, axial ligands interact according to the "extended first solvation shell" paradigm.This notion is further supported by the asymmetry in the average distance of the two axial MeCN ligands with respect to Pd 2+ (Fig. S9 and S10).Additionally, a detailed analysis of the axial coordination pattern (see section §S3.1) suggests that the preferred average structure of the complex is not octahedral but square pyramidal, with the ratio between the latter and the former being 7:4.

Free energy barrier of ligand exchange
The ability of the MACE potentials to describe ligand exchange around the metal centre relies on accurately describing the dif-ferent coordination states and their exchange mechanisms, ideally leading to accurate kinetics.The latter has been difficult to achieve with classical force fields. 59,60 discussed previously, ligand exchange in Mg 2+ complexes is suggested to follow a dissociative or interchange-dissociative mechanism, which proceeds through an intermediate or transient structure with a lower coordination number.Solvent dissociation from the first solvation shell of the [Mg(H 2 O) 6 ] 2+ thus represents the rate-limiting step of the process.The free energy barrier associated with the dissociation of one solvent molecule was obtained from the potential of mean force (PMF) using US with 48 windows (Fig. 4a).The PMF shows a minimum at 2.1 Å, in agreement with the value obtained from the RDF (Fig. 3a).A second shallow minimum is located around 4.25 Å, corresponding to the position of the second solvation shell, which is approx.2 kcal/mol higher than the first minimum.This indicates that the sampling of the region where water leaves the first solvation shell is not fully converged.However, the two repetitions of the US confirm that the free energy barrier is not affected (see Fig. S13).The representative structures of both minima and the associated transition region are depicted in Fig. 4. Analysis of the trajectories confirms that the MACE potential correctly restores the octahedral geometry of the complex.The barrier of pulling the water molecule away from the first solvation shell is 7.6 ± 0.15 kcal/mol, with the peak located at a distance around 3.1 Å.The predicted barrier is 1.9 kcal/mol lower than the experimental value, 9.5 kcal/mol, 4,11 implying a higher rate of water exchange around the Mg 2+ .
40,41 To determine the coordination number that corresponds to the TS, the fluctuations of the coordination number and the Pd−N bond lengths were analysed in the umbrella sampling trajectory that corresponds to the TS region in the PMF (Fig. S16).The coordination number fluctuates around a mean value of 4.96, which yields a barrier of 15.7 kcal/mol in the PMF (black dot in Fig. 4b).This value is in very good agreement with the experimental values from two independent NMR The reaction mechanism of the ligand exchange is illustrated in detail by the representative snapshots in Fig. 4b.The solvent association proceeds via the formation of a square pyramid, where the axial Pd−N bond shortens as the system approaches the TS region, while one of the equatorial Pd−N bonds progressively elongates.This process leads to the formation of a trigonal bipyramidal TS, with an equal length of two Pd−N bonds involved in the ligand exchange.Taken together, the PMF from the US simulations confirms the associative nature of the ligand exchange mechanism.Furthermore, the activation barrier is in excellent agreement with the experimental values, supporting the notion that the MACE potential fitted to the hybrid DFT reference accurately describes the PES of the system and allows for a realistic description of the dynamics of the ligand exchange process.

Conclusions
In this work, we present computational strategies to build training data sets for modelling ligand exchange processes of divalent metal cations in explicit solvents with MLPs.Using Mg 2+ in aqueous solution and Pd 2+ in MeCN as model systems and illustrative examples, we demonstrate the capability of the MACE potentials to reproduce the total energies and forces of the solvated metal cations.Furthermore, the MLPs trained on cluster data can be used in the condensed phase simulations with periodic boundary conditions.The MACE potentials yield metal ion-solvent RDFs in excellent agreement with experimental data, confirming the capacity of MACE to capture the structure of the polarised solvent shells around the cations.Moreover, we demonstrate the ability of the MACE potentials to model changes in the coordination shells of the metal cations, allowing for a structurally and energetically realistic description of different ligand exchange mechanisms in complex liquid environments.More generally, we show that the active learning strategy combined with MACE potentials allows the generation of accurate and data-efficient MLPs that are suitable to model changes in the coordination chemistry of charged species in solution.While further work is needed to automate the selection of accurate parameters suitable across different metals, this study provides a robust computational framework for preparing data-efficient models that accurately describe metal-ligand interactions, paving the way to modelling increasingly complex systems, such as metallocages and catalysts.

Production MD
Mg 2+ cation in water was simulated in a periodic cubic box of 16.3 Å containing 1 cation and 145 water molecules.The initial structure for the dynamics was obtained by classical dynamics (see SI §2.3.1).The Pd 2+ system was simulated in a periodic cubic box of 24.0 Å containing 1 metal atom and 159 acetonitrile molecules.The starting configuration of the system was built in three steps.Firstly, an NPT-equilibrated box of acetonitrile was obtained with the force field parameters from Caleman et al. 106 .Secondly, the xTB-optimized structure of the [Pd(MeCN) 4 ] 2+ complex was solvated with the acetonitrile box.Lastly, the whole system was energy-minimized with Gromacs (v2019.1). 107The radial distribution functions for both the Mg 2+ and Pd 2+ systems were computed from 500 ps NVT simulations with the MACE potentials, with the first 50 ps used for equilibration and skipped from the analysis.The equations of motion were integrated by the i-PI driver, with MACE potential evaluated by the MACE-Atomic Simulation Environment (ASE) calculator using ASE v3.23.0b1. 109All MD simulations were propagated with an integration time step of 0.5 fs.MD in the NVT ensemble was thermostatted at 300 K by a stochastic velocity rescaling thermostat with a coupling time constant of 100 fs. 110

Free energy computations
The free energy profiles of the metal-ligand exchange reactions were evaluated by umbrella sampling (US) simulations using the i-PI driver combined with the Plumed v2.9.0 library. 111,112The PMFs were constructed using the Weighted Histogram Analysis Method (WHAM) code v2.0.11. 113r Mg 2+ in water, the Mg•••O distance was chosen as a recoordinate (RC).The starting structures for each window were generated by a steered MD, pulling a water molecule from the first solvation shell to a distance ranging from 1.5 to 7 Å.The US covered the distance from 1.5 to 6.0 Å, split into 48 windows with a spacing of 0.075 Å.In each window, the trajectory was propagated for 1 ns, with the first 50 ps used for equilibration and skipped from the analysis, corresponding to 45.6 ns of sampling time.The US windows were propagated independently at 300 K.A harmonic umbrella restraint with force constant 500 kcal/(mol•Å 2 ) was employed in the windows from 1.5 to 3.00 Å, while the force constant was lowered to 250 kcal/(mol•Å 2 ) in the region from 3.00 to 6.00 Å.The PMF was computed as an average from 2 US runs with random seeds for generating the initial velocity distributions, with an uncertainty estimated as a standard deviation.The final PMF was corrected by the entropy term 2k B T ln(r/r 0 ) that accounts for the increasing volume of configuration space with increasing distance r. 114,115 For Pd 2+ in acetonitrile, the coordination number was chosen as a collective variable for the US (further details in section SI §S4.2).The starting structures for the US runs were generated by steered MD, during which one of the two MeCN axial ligands was pushed towards the Pd centre, guiding the system towards the ligand exchange event through the formation of the pentacoordinated TS.In the US, the coordination number was varied

Fig. 1
Fig. 1 Training data and active learning workflow: a) Training subsets for Mg 2+ in aqueous solution, b) training subset for Pd 2+ in acetonitrile (MeCN), c) Scheme of the active learning workflow used to train the machine-learning potentials (MLPs).

OpenFig. 3
Fig. 3 Simulation boxes, radial distribution functions g(r) and coordination numbers N(r) of the metal complexes in solution.a) Mg 2+ in aqueous solution, b) Pd 2+ in MeCN.

FaradayFig. 4
Fig. 4 Potential of mean force (PMF) profiles of the two ligand exchange processes for a) [Mg(H 2 O) 6 ] 2+ (the solvent molecule exchanged is coloured blue) and (b) [Pd(MeCN) 4 ] 2+ , where the black dot indicates the energy at the TS.