Ion Association in Lanthanide Chloride Solutions

Abstract A better understanding of the solution chemistry of the lanthanide (Ln) salts in water would have wide ranging implications in materials processing, waste management, element tracing, medicine and many more fields. This is particularly true for minerals processing, given governmental concerns about lanthanide security of supply and the drive to identify environmentally sustainable processing routes. Despite much effort, even in simple systems, the mechanisms and thermodynamics of LnIII association with small anions remain unclear. In the present study, molecular dynamics (MD), using a newly developed force field, provide new insights into LnCl3(aq) solutions. The force field accurately reproduces the structure and dynamics of Nd3+, Gd3+ and Er3+ in water when compared to calculations using density functional theory (DFT). Adaptive‐bias MD simulations show that the mechanisms for ion pairing change from dissociative to associative exchange depending upon cation size. Thermodynamics of association reveal that whereas ion pairing is favourable, the equilibrium distribution of species at low concentration is dominated by weakly bound solvent‐shared and solvent‐separated ion pairs, rather than contact ion pairs, reconciling a number of contrasting observations of LnIII–Cl association in the literature. In addition, we show that the thermodynamic stabilities of a range of inner sphere and outer sphere LnClx(3-x)+ coordination complexes are comparable and that the kinetics of anion binding to cations may control solution speciation distributions beyond ion pairs. The techniques adopted in this work provide a framework with which to investigate more complex solution chemistries of cations in water.


Introduction
Lanthanides (Ln), which form the majority of the rare-earth elements,a re used extensively in technological devices due to their magnetic, electronic and opticalp roperties. [1] Coordination complexes of the lanthanides are particularly important in display technologies [2] anda sc ontrasta gents in medical imaging. [3] Furthermore, there is growing interest in lanthanide complexes as agents in the treatment of cancera nd in biology. [4,5] Given that the security of supply of the lanthanides is precarious, [6,7] significant effortsa re being made to identify methods to efficiently extract, processa nd recycle them. [8] Separating out different lanthanides is ap articularc hallenge that is usually performed using solvente xtractionm ethods that tend to be expensive, hazardous and environmentally damaging. [9] Solvent-based processing is also used extensively in nuclear-waste management for lanthanide and actinide separation. [10][11][12] For all of the above reasons, ar obust understanding of the solution chemistry of the lanthanides is of great importance.
In aqueous solutions, the lanthanides are principally found in the Ln III oxidation state. Ln III ionic radii decrease on moving across the series due to increasing nuclear charge, and this re-sults in ac oncomitantd ecrease in the number of water molecules in their aquo complexes in solution, affecting their transport properties. [13] It is widely accepted that the most likely water coordination number (CN) in the first shell for the light Ln III is nine, with complexes adoptingatricapped trigonalprism (TTP) geometry. [14] The CN decreases to eight for heavy Ln III ,f or whichs quare-antiprism (SAP) complexes are usually observed. Althoughe arlier experimentsr eported water CNs for Ln III in the middle of the series as eight [15,16] or nine, [17] it is now widely accepted that the average CNs fall between these values. [18,19] Water-exchange rates between the first and second Ln III coordination spheres appear to increase across the series and reach am aximum aroundg adolinium before subsequently decreasingf or the heavier lanthanides; [14] however,e xperimental data for water exchange is incomplete for all Ln III .T he mechanism for water exchange was proposed to be associative-activated exchange for eight-coordinatec omplexes and dissociative-activated exchange for nine-coordinate complexes. [14,20] Given that the relaxivity of water surrounding Gd III is extremely important for its use as ac ontrasta gent in magnetic resonance imaging, there is particulari nterest in the water-exchange rates and protonr elaxation for solution complexes of this lanthanide. [3,21] Exchange of the water moleculess urrounding Ln III with anions may occur in common electrolyte solutions to form ion pairs, trimers, etc. Of particulari nterest is the formation of [LnCl x (OH 2 ) y ] (3Àx)+ + ,not simplyb ecause Cl À is ac ommona nion in the laboratory,b ut also duet ot he high Cl À concentrations, c(Cl À ), that are found in many naturald eposits containing elevated levels of Ln III . [22,23] Stability constants, b,p rovide the activitieso fa ssociated species relative to those of individual constituents; for example, for Ln III chloride coordination complexes, b x Cl (where x is an integer > 0) indicate the dominant speciesassociated with the following equilibria, For the formation of an ion trimer (x = 2), where a refers to the activities of species and b refers to ac umulative constant which may comprise multiple stepwise-equilibrium constants.
In the formation of the ion pair in which the ions are in direct contact (the contact ion pair (CIP)) from dissociated ions in solution, the system must first sample the solvent-separated ion pair (SSIP) and the solvent-shared ion pair (SShIP) states, in which one or two shared water shells separate the ions, respectively,a ss ummarised in Figure 1. If K is the equilibrium constant for the stepwise formation of each type of ion pair (e.g.,f rom the left to the right of the Scheme in Figure 1), then in theory, b 1 Cl ¼ K SSIP ð1 þ K SShIP ð1 þ K CIP ÞÞ,b ut in practice, this is only true if the energetic barriers between states can be overcome on the timescales of the experiment (i.e.,t he system achievesatrue equilibrium).
In as tudy by Gammons et al., [24] the stabilities and concentrationso fN dCl x speciesw ere determined indirectly by considering the solubility of AgCl in solutionsw ith varying concentrations of Nd 3+ + and Cl À above 40 8C. Extrapolating from the formation constants at higher temperatures,t hey determined that log 10 (b 1 Cl )a tT = 25 8Cw as 0.06 AE 0.5 and that no further association by Cl À is found at this temperature. Both b 1 Cl and b 2 Cl increased by several orders of magnitude above 50 8C. Also, by indirect measurement of solution species, Luo andB yrne [25] found little variation in b 1 Cl across the lanthanide series. The mean log 10 (b 1 Cl )w as 0.65 AE 0.05 at T = 25 8Cb ye xtrapolating to infinite dilution.
In line with the conclusions from earlier work by Mundy and Spedding, [26] Allen et al., [27] using extended X-ray absorption fine-structure (EXAFS) spectroscopy,d etermined that outersphere Ln III chloride complexesa re formed at low c(Cl À ). When c(Cl À )e xceededa bout1 0mol dm À3 ,i nner-sphere complexes with average coordination 2.1 (La) to 1.1 (Eu) wereo bserved. The decreasing coordinationa tc onstant ionic strength, I,s uggests that the stability of chloride complexes decreases across the series. This is also consistentwiththe earlier data of Martell and Smith [28] in which log 10 (b 1 Cl )w as reported to be 0.48-0.23 for La-Lu (where I = 0.1mol dm À3 and T = 25 8C). The Lawrence Livermore National Laboratory thermodynamic database [29] provides log 10 (b 1 Cl ) = 0.3086, log 10 (b 2 Cl ) = 0.0308 and log 10 (b 3 Cl ) = À0.3203f or the association of free chlorides to Nd.
Given the discrepancies between different experiments,i ti s reasonable to assumet hat ion pairing is controlled by kinetic factors, potentially linked to water removal aroundt he cation. It is difficult to directly comparet he relative stabilities of different LnCl x species in experiments due to the changes in water activities associated with what are often large changes in I. Furthermore, the methods used to estimate formation constants are inconsistent across the experimental set. Simulations at the atomicl evel offer an invaluablea lternative approach and have made significant advances in our understanding of Ln 3+ + (aq) over the previoust hree decades. We note though that even here there are significant challenges associated with comparing experimental measurements, such as equilibrium constants, to calculated values from simulations. An umber of computational studies of ion pairing have made progressi n dealingwiththis challenge. [33][34][35] Early classical molecular dynamics (MD) simulations by Meier et al. [36] showedt hat CN in the first La III coordination shell was concentration dependent( waterC Na nd chloride CN decrease and increase, respectively,w ith increasing c(Cl À )), but their values of averagec oordination are now known to be too large. Kowall et al. [20,37] improved on this by including some polarisability in their water model which was suggestedt ob ei mportant to accurately capturew ater-exchange dynamics. Floris and Ta ni [38] developedp otentials based on ab initio calculations and found the capped square antiprism (CSQA) as an alternative structure to aT TP for Ln III nona-aquo complexes.S ignificant understanding of Ln 3+ + (aq) structurea nd dynamics has been provided by Duvail et al.,p articularly when simulations were combined with spectroscopic experiments. [18,39] They included explicit polarisability of water which was claimedt ob e crucial for the accurate prediction of the structure and dynamics of the first Ln coordination sphere. [40,41] Furthermore, simulationsw ith chloride, perchlorate and nitrate allowed the authors to quantify the affinity for anion association to Ln III by potentialofm ean force (PMF) calculations. [42] Villa et al. [43] developed af lexible, polarisable force field for Ln III in water based on ab initio calculations. Reasonable agreement with experiments was found in the structure and dynamics of Ln III aquo complexes. To achieve water-exchange rates comparable with experiment,h owever,r equired reducing the polarisability from that calculated using ab initio measurements. Modelling explicit polarisability is thought to be neces- Steps in the formation of ac ontact ion pair( CIP) from solvated "free" ions in solution. As the cationapproaches the anion, solvent molecules surrounding ions are displaced to form the solvent separated ion pair (SSIP) and the solvent-shared ion pair (SShIP) beforethe CIP in as eries of stepwise equilibria. Generic cations, anions and solvent molecules are shown by the green,o rangea nd blue circles,r espectively. sary to accurately capture ion solvation energies. [44,45] Te rrier et al., [46] using density functional theory (DFT), showed that the polarisation of water in the first La III shell led to a + 0.5 Ds hift in the dipole momento fw ater molecules compared with those in the bulk. Ar ecent ab initio study [47] shows that this polarisation is dominated by charge-dipole interactions for metal ions in water and that isotropic polarisability is capable of accurately modelling these systems. Migliorati et al. [48,49] have recently published pairwise intermolecular-potentialp arametersf or all Ln III with water.B ased on X-ray spectroscopy data, with an on-polarisable force field and rigid water,t he modelsw ere able to accurately reproducet he structuralp roperties of Ln III aquo complexes. Qiao et al. [50] also developed a non-polarisable force fieldb ased on CHARMM and compared it to the polarisable force field of Marjolin et al. [45] The solvation energiesf or the non-polarisable model werefound to be as accurate as those calculated with explicit polarisation included, suggesting that polarisability is not crucial forag ood quality model.
Despite the previous work investigating Ln III aqueous solutions, inconsistencies remaini no ur understanding of these systems, particularly with regard to cation speciation. The mechanisms by which ions associate in solutionr emainu ncertain and the thermodynamic stabilities of coordination complexes are unclear.I nt he present study,M Ds imulations have been adopted to explore the structurea nd stabilityo f [LnCl x (OH 2 ) y ] (3Àx)+ + in aqueouss olutions. Nd, Gd and Er have been simulated in chloride solutionso ver ar ange of concentrationsu sing both density functionalt heory (DFT-MD) and classical empirical potentialM D( C-MD) methods. These calculations show that the effect of ion size controlsb oth the structure and association mechanisms of coordination complexes in solution.

Methods
We have parameterised an umber of interatomicp otentialenergy functions to simulate Ln III in aqueous solutions. The force-field detailsa re provided in Ta ble 1. Both DFT-MD and C-MD have been adopted in this study,w ith the DFT-MD simulations and experimental data being used to benchmark the C-MD calculations.

System Preparation
Simulations werep erformed for one Ln III in 5550, 555 and 111 water molecules (b(Ln) = 0.01, 0.1 and 0.5 mol kg À1 ). The number of Cl À ions was varied from zero to three depending on the particular study.
Unless otherwises tated, simulations were performed at a constantm ass density and 298 Kw ith cubic cells and 3D periodicity.T he mass densitieso fs olutionsw eres et according to the empirical equation defined by Speddinge tal., [51] In the above equation, d and m are the target mass density and concentration of the solution, whereas d 0 is the mass density of water and A i are as eries of density parameters fitted to experimental datat hat can be found in Ref. [51] and in Ta ble S1 in the SupportingI nformation. We checked the validity of these densities in simulations at constantt emperature and pressure.
To calculate the structural properties of Ln III complexes in water,M Dw as performed for approximately 14-17.5 ps and 5nsatb(Ln 3+ + ) = 0.5 and0.01 mol kg À1 in DFT-MD andC-MD, respectively.S imulations were performedw ith three additional Cl À ions in aqueous solutiona tb(Cl À ) = 1.5 (DFT-MD) and 0.03 mol kg À1 (C-MD) as well as in pure water.W ater molecules in the first Ln III coordination shells were assigned using ad istance criterion, in which the cutoff for LnÀOw as 3.5 ,a si nformed by respective time-averaged radial-distribution functions (RDFs). We chose to initiate C-MD simulations with either eight (Gd and Er) or nine (Nd) water molecules in the first shell and the simulations werea llowed to equilibrate over 2ns. In the case of the DFT-MD calculations, the first Ln 3+ + shell was set to both eight and nine water molecules in different simulations for Nd and Gd, (toe nsure our sampling was not restricted by the short trajectories of DFT-MD as the water-exchange time is of the order of the total DFT-MD simulation time) [14] and eight water molecules for Er. Dipole moments from DFT-MD were calculated using the maximally localised Wannier function( MLWF) formalism [52] which providesapictureo ft he electron distribution around atoms. [53] Thisf ormalism was appliedt oc onfigurations that were extracted every 10 fs from DFT-MD simulations. Thus, a minimum of 100 configurations were analysed for each study of Ln III complexes in water.
Umbrella Sampling (US) [54] was appliedt om easure the freeenergy change for the reversible binding of Cl À with Ln 3+ + (that is, Eq. (1) with x = 1) using the Plumed plugin [55] in C-MD simulations at b(Ln 3+ + ) = 0.01 mol kg À1 .I no rder to achieve ah igh level of accuracyf or the energy barriers between states andt o gain mechanistic insight into ion association, we sampled a two-dimensionalr eaction coordinate as defined by two collective variables (CVs). Harmonic-potential energy biases were imposed to restrain the distance between Ln 3+ + and one Cl À ion in the range 2.5-16 at either 0.25 or 0.5 intervals. Twoa dditional chloride ions were restrained beyond two water shells away from Ln 3+ + (approx.1 0a nd 15 separation).I na ddition, the water oxygen (Ow) coordination number in the first Ln shell, as defined by the continuous function, was restrained between values of 6.5 and 9.5 at intervals of approximately 0.5. In Equation (4), N is the number of water oxygens in the simulation and r i are the LnÀOd istances.T he remaining parameters are set such that the functions moothly goes from one to zero within the bounds of the first and second Ln solvation shells (for more details see the Supporting Information, Figure S1 and Ta ble S2). Windows were added or removedd epending on the system under investigation;p robability distributions in the 2D reactionc oordinate were monitored to ensure good overlap between CV distributions in adjacent US windows in the sampled space. See Ta [56] was subsequently used to generatep otentialofm eanforce (PMF), W,e nergy maps. Metadynamics [57] calculations were performed to measure freeenergy differences between [LnCl x ] (3Àx) + species where x = 0-3 and b(Ln 3+ + ) = 0.1 mol kg À1 using C-MD. At hree-dimensional reaction coordinate was adopted, in which the three sampled CVs were S LnÀOw for O( of water) and S LnÀCl1 for Cl coordination number in the first Ln shell, respectively,a nd also S LnÀCl2 for Cl coordination number in the first and second Ln shells combined. Coordinationo fL nb yc hloride was defined as in Equation (4) and the parameters for each CV can be found in the Supporting Information in Ta ble S2. Three metadynamics calculations for each Ln system were initiated from different starting configurations, one for each of x = 0, 1a nd 2. Gaussian-shaped bias potentials were deposited every 1000 steps with ah eight of 0.5 kJ mol À1 .T he width of Gaussians in the CVs forc hloride coordination was 0.015, whereas for Ln-water coordination, this was 0.033.S imulations were performed for approximately 100 ns and the resulting free-energy surfaces were generated by summing all of the deposited time-dependent bias. Changesi nf ree energiesw ere then calculated as the average difference in the free energies between the same regions of CV space from three independent calculations with associated standard deviations.

DFT-MD computational details
Simulations were performed using the ViennaA bi nitio Simulation Package (VASP), [58,59] employingt he projector augmentedwave method. [60,61] Pseudo-potentials were generated using valence configurations of 5s 2 5p 6 5d 10 6s 2 for Nd, Gd and Er (with f electrons kept frozen in the core), 3s 2 3p 5 for Cl, 2s 2 2p 4 for O and 1s 1 for H. Scalar relativistic effects were accountedf or,b ut spin-orbit coupling was neglected. Simulations were run at T = 341 Ku sing the optB88 [62,63] van der Waals density functional. [64][65][66] Previous studies have identified this temperature as the most appropriate for DFT-MD simulations of liquid water (the enhanced temperature compensates for over-structuring of the water by the functional and gives ag ood description of the liquid water structure and self-diffusion coefficient at ambient conditions) [67] and is the best choice for simulations of large systems. [68,69] The kinetic-energy cut-off for the planewave expansion was 600 eV andBrillouinzone sampling was restricted to the G-point. The break condition for the electronic self-consistent loop was set to 10 À5 eV.M olecular dynamics simulations were performed using at ime-step of 0.5 fs, making use of the Nóse-Hoover thermostat. [70,71] The atomic mass of hydrogen was set to 2amu. whicha llows us to use a longer time step.

C-MD computationald etails
The SPC/Fw flexible three-point water model [72] was adopted because this has been shown to perform well in simulationso f metal ions in water [73,74] as well as accurately capturing liquidwater self-diffusion coefficients and dielectric properties. The mass density of SPC/Fwa t2 98 Ka nd 1atm is 1012 kg m À3 , therefore d 0 in Equation (3) was set to this value. The Ln 3+ + ÀO (water) intermolecular potentialp arameters were calculated by first optimising crystal structures of Nd 2 O 3 , [75] Gd 2 O 3 [76] and Er 2 O 3 [77] in the GULP simulation package [78] using the force field developed by Lewis and Catlow. [79] Given that no Er 3+ + ÀO 2À interatomic potentiali si ncluded in the originalf orce field, we generated Buckingham potential-energy A and 1 parameters by interpolationf rom the parameters available for other Ln 3+ + À O 2À interactions. The methodo fS chrçder et al. [80] was then used to fit the interaction potentialb etween Ln and O. This involveds caling the Ln and Oc harges (q Ln = 2.013 e, q O = À1.354 e), accordingt ot he chargeo nOin the SPC/Fw model, before optimising the Buckingham potential parameters (see Ta ble 1) to reproduce the crystal structure with thesen ew charges. The procedure ensures that the Pauli repulsion modelled between atoms in the solution phase is at least consistent with that for the crystalline phase. To model the Ln-Cl intermolecular interactions, NdCl 3 ,G dCl 3 and ErCl 3 crystal structures from experiments [81,82] were used to fit Buckingham A and 1 parameters while retaining the crystal symmetry.C rystal geometries were allowed to relax duringt he fitting except in the case of ErCl 3 for which geometriesw ere fixed because relaxation led to changes in the Er-Cl coordination environment (although the total coordination number remained constant).I ti sk nown that ErCl 3 is unstable due to its highly hygroscopic nature, and so there is limited knowledge on the most stable crystal structure. LnÀCl distances were monitored to ensure that the mean first coordination-shell distances in the final structure were within 0.1 of the starting structure.U sing this approach, it is possible to generate a spectrum of A and 1 pairs. Linear fitting of 1 versus A provides af unction with which to calculate the intermolecular energies for all A and 1 parameters that conserve the first shell LnÀCl coordination distances. We chose A and 1 based on the maximum potentiale nergies in the short-range Buckingham potential at the mean LnÀCl distance following anumber of tests. Althought his choice does provide an upper bound to the interaction energies, changes in Buckingham energies were within thermale nergy at 298 Kw hen comparing the upper and lower bounds of the identified parameter set from crystal-structure fitting.
Simulations were performed using the DL_POLYC lassic package. [83] Atom trajectories were obtained using aV erlet Leapfrog algorithm with a0 .5 fs timestep. Af ive-chain NosØ-Hoovert hermostat [84] with a0 .1 ps relaxation time was employed to maintain the target temperature at ac onstant volume. US calculations were performed with one Ln 3+ + ,t hree Cl À and 5550 H 2 O( i.e., b(Ln 3+ + ) = 0.01 mol kg À1 ), whereas metadynamics werec onducted with one Ln 3+ + ,t hree Cl À and 555 H 2 O( i.e., b(Ln 3+ + ) = 0.1 mol kg À1 )a nd the cubic-simulation cell parameters were set accordingt ot he target density [see Eq. (3)].T om easures olvation enthalpies,s imulations containing one Ln 3+ + and 5550 H 2 Oa tc onstant temperature and pressure were performed. Here, temperaturea nd pressure were constrained using aN osØ-Hoover thermostat and barostat with 0.1 ps and1 .0 ps relaxation times, respectively.S hortrange intermolecular interactions were truncateda tadistance of 9 and as mooth particlem esh Ewald [85] algorithm with 10 À7 precision was used to calculate electrostatic potential energies. In systems with an et charge, an eutralising uniform background chargew as appliedu sing the correction of Fuchs. [86] Results and Discussion Ln III aquo-complex properties Here we describe the resultsf rom C-MD and DFT-M Ds imulations with Ln 3+ + and 3Cl À ions in aqueous solution.F igure 2 highlightst he very good agreement in the LnÀOR DFs between the two types of simulation. The DFT-MD RDF peaks tend to be narrower than the corresponding C-MD ones, which can be explained by stronger Ln III ÀOi nteractions and/or by the shorter duration of the simulation runs. Indeed, water-molecule exchanges between the first and second Ln III coordination shells are observed less frequently in DFT-MD calculations due to the timescale of these processes being at least on the order of the simulation times. [14] The positions of the maxima in the first and second peaks of the LnÀOR DFs (d LnÀO1 and d LnÀO2 )i nF igure 2a re given in Ta ble 2. The parameter d LnÀO1 increased by approximately 10 % on moving from the heaviest to the lightest lanthanide. The RDF peaks for second-shell water were shifted to slightly larger distances for Gd and Er.P eak positions and coordination numbers generally comparew ell with other theoreticala nd experimental studies (see Ta ble 2). Consistent with earlier studies [87,88] we see little evidence for athird solvent shell.
No inner-sphere Ln III chloride complexes formed in solution. RDFs (see the Supporting Information, Figure S3) show that in all cases, the shortest Ln-Cl distances were beyond4.W hen chloridew as restrained to the innermost Ln III coordination sphere, the cation-anion distance ( Figure S1) was consistent betweent he DFT-MD and C-MD. In the second shell,C l À was most likelyt ob efound in the hydrogen coordination sphere surrounding Ln III .V ery little change was seen in Ln-water RDFs or CN in the second Ln III shell when comparedt os imulations in pure water.T he largest difference between the DFT-MD and C-MD is in the coordination distances of water hydrogen to Ln III in the first shell. For example, in NdCl 3 (aq) simulations, the maximum in thef irst peak of NdÀHR DFs was at 3.13 and 3.24 in DFT-MD and C-MD, respectively.T his is most likely due to the explicit polarisabilityi nD FT leading to stronger water-Ln interactions and ad ifferent orientation of water moleculesi nt he DFT-MD compared with the C-MD, as discussed in further detail later in the text. By the second coordination Average values for the water CN in the first and second Ln III coordination shells (CN 1 and CN 2 )a re provided in Ta ble 2. CN 1 varies linearly when modelled with the classical empirical potential. For the largest cation,t here is close agreementb etween CN 1 in all of the presented studies. The DFT-MD CN 1 for Gd appears to be too low when compared with C-MD and other studies. Given the relatively long timescales for water exchange around Ln, DFT-MD is unlikely to sample ar epresentative range of equilibrium states for Ln coordination complexes in water.T he value of CN 1 for Er is not clear from the wider literature, and both the DFT-MD and C-MD values fall within the range reported. Figure 2s hows the probabilitiesf or the number of water molecules in the first shell when the truncation distance for oxygen coordination was 3.5 .T here is a2% probability of octa-aquoc omplexes in the case of Nd and the coordination is dominated by nine water molecules. For Gd and Er,t he octa-and nona-aquo complexes are both observed with quite high probabilities. Largerd eviations are found in both d LnÀO2 and CN 2 in this work and the wider literature. Given that the structuring of water in Ln III shells is greatlyr educed beyondt he first shell, it is likely that solution conditions will affect the average number of water molecules even at these distances. Figure S4 in the Supporting Information shows O-Ln-O angles and Ln-O-M angle probabilityd istributions in which M is the bisectorb etween the two water OÀHb onds;h ence, the Ln-O-M angle quantifies the relative tilt of water molecules surrounding the cation. Bimodal distributionsa re observed in both DFT-MD and C-MD for Nd, consistent with earlier work. [41,50,88,89] As discussed by Qiao et al., [50] the peaks centred around7 0 8 and 1358 represente ither the TTP or the gyroelongated square-antiprism (GySQAP) geometries with slight deformations. The TTP geometries dominated our simulations of nona-aquo Nd complexes but we did find as maller population of GySQAP complexes. Representations of these structures are shown in Figure 3.
With decreasing lanthanide-ions ize, O-Ln-O peaks shift closer to the positions of peaks representing aS QAP geometry (758 and 1428). The shift was more pronounced in DFT-MD simulationsi nt he case of Gd, because the starting point here was an octa-aquo complex with SQAP geometry.I na ddition to these shifts, ap eak centreda round1 128 was found for the heavierL n III in other simulations. [41,50] This was observed in our studies as ab roadening of the second peak towards smaller anglesa nd was most notable in C-MD simulations of Er 3+ + .
The Ln-O-M angle distributions in Figure S4 (Supporting Information)s how ah igh probability for angles close to 1708 for all simulated lanthanidesi nC -MD simulations. This indicates that the majority of water molecules in the first coordination sphereb ind to Ln III in such aw ay to reduce any electrostatic repulsion with positively charged water-hydrogen atoms. The probability was reduced in the case of DFT-MD simulations. Althought he maximum probability in Ln-O-M angles was again similarf or all Ln and close to 1708,p eaks were much broader with angles of~1408 showings ignificant probability.Asimilar distribution was previously calculated by Te rrier et al. [46] for Table 2. The mains tructural,d ynamical and energetic features of Ln III in this worka nd other experimental and modelling studies.L n-O water distances in the first (d LnÀO1 )a nd second( d LnÀO2 )L ncoordination shells are given in units of .A veragec oordination-number valuesw ere calculated by counting the number of water oxygensi nt he first (CN 1 )a nd second (CN 2 )h ydration spheres whereadistance criterion was used-informed by the positions of minimai nR DFs.W ater-exchange rate constantsb etweent he first and secondL nc oordination sphere, k ex ,a re given in ns À1 units. H solv are in units of kJ mol À1 (see the discussion in the text regarding the recalculation of data from Ref. [95]).  [18] 2.53 9.0 MD [41] 2.48 9.0 4.63 19.2 0.67 MD [48] 2.52-2.538 .8-8.9 XRS [101] 2.51 8.9 MD [43] 2.63 8.9 0.55-0.83 À3524 17 ONMR [14] ! 0.5 MD [44] À3429 Expt.r ecal. [95] À3501 Gd III DFT-MD [ [18] 2.46 9.0 MD [41] 2.39-2.448 .7-9.0 4.55-4.611 8.9-19.22 .35-3.94 MD [48] 2.45-2.468 .6-8.8 17 ONMR [14] 0.83 MD [102] 2.44 8.6 4.65 18.1 2.69 MD [43] 2.55 8.4 1.61-1.81 À3659 MD [44] À3617 Expt.r ecal. [95] À3601 Er III DFT-MD [ [48] 2.39 8.0-8.3 XRS [103] 2.37 8.2 17 ONMR [14] 0.13 MD [44] À3740 Expt.r ecal. [95] À3726 [a] This work. La 3+ + using Car-Parrinello MD simulations. Given that the first peak in Ln-H RDFs (see Figure S3 in the SupportingI nformation) was at shorter distances in DFT-MD, it is sensible to conclude that water molecules in the first shell show someh ydrogen-bonding capacity with each other as wella sw ith water molecules in the second shell, explainingt he smaller value of CN 2 by 10-20 %i nD FT-MD than in C-MD. Given the above result,w ec ompared the dipole momento f water molecules, m H 2 O ,i nt he first Ln III coordinations pheres to those of all other water molecules beyond the second shell, that is, water in bulk solution.F igure S5 (Supporting Information) provides dipole-moment probability distributions.T he dipole momentf or bulk SPC/Fw water is 2.4 Dw ith as tandard deviation of 0.16 D, smaller than that measured for bulk liquid water in experiments (2.9 AE 0.6 D) [90] and the value estimated by Car-Parrinello MD (3.0 D). [91] Classicall iquid-water models often underestimate experimentally determined m H 2 O measurements. [92,93] m H 2 O in DFT-MD for bulk water was 3.1 AE 0.3 D.
m H 2 O distributions for water molecules in the first Ln III coordination shells in C-MD were shiftedt og reater values than those found for bulk water (5-7 %f rom Nd to Er). The biggest increase was for the lanthanide with the highest charged ensity. The distributions in the second solvation spheres were close to those for bulk water.T he trend in first-shell water perturbation was also observed in DFT-MD simulations. Here, m H 2 O distributions were more perturbed than for the C-MD calculations. Indeed,i nt he case of Nd, Gd and Er,s hifts to greater mean values by approx. 0.37, 0.54 and 0.52 D, respectively,w ere found to be in agreement with the results from DFT simulations of La 3+ + (aq). [46] Water-molecule exchange between the first and second Ln III coordination spheres was monitored throughout the C-MD simulations. It was not possible to measure this with anyr easonable accuracy from DFT-MD, given the slow dynamics of the process. Any water-oxygen atomst hat were within 3.5 of Ln III were considered to be inside the first shell. By measuring the time that water molecules reside in the first shell after first entering( t r )w eg enerated water-residence probability distributions, P(t r ). Exponentials of the form, P(t r ) = p(t 0 )exp(Àt/t r ), were fitted to the data to estimate water mean residence times (t r ). The exchange rate (k ex )i so btained as the reciprocal of t r and these values are provided in Ta ble 2. Given the inherentbias associated with defining coordination according to af inite distance criterion, we also measured residence times in which the truncation distance for oxygen coordination to Ln III was the maximum in the second peak of LnÀOR DFs (see Figure 2) to provide aliberal interpretation of first-shellcoordination. Only water-molecule exchange eventst hat persisted for at least 10 ps were considered in all our analyses.
Exchange rates forN d III compare reasonably well to those from other theoretical studies and from 17 On uclear magnetic resonance (NMR) experiments; however,f or Gd III and Er III , k ex values are considerably larger than in experiments. It is worth noting the wide variation that is found in all of the calculated k ex values in Ta ble 2a nd that MD data from others tudies also predict more frequentw ater exchange between the first and second Ln III shells than experiments.A sr eported in Helm and Merbach, [14] the general consensus from 17 ONMR studies is that k ex increases from La, reaching am aximum for lanthanides in the middle of the series before subsequently decreasing. This is explained by the transition from predominantly nine to eight-fold coordination that is observed in the middle of the series.W eo bserve only increasing exchanger ates across the lanthanides and so it may be either that simple models cannot capturet he kinetics of water exchange fully or that we are seeing ar eflection of uncertaintya round the experimental values. Duvaile tal. found that the mean residence time for water molecules surrounding Gd decreased relative to that of Nd, but for Er this was approximately the same as Gd. [41] The number of exchange events( water molecules either leaving or enteringt he first solvation sphere) per nanosecond further reflects the lack of am aximum water-exchanger ate in the middle of the lanthanide series:a pproximately 20, 90 and 106 ns À1 ,for Nd, Gd and Er,r espectively.
The diffusion coefficients, D,f or Nd, Gd andE ri nS PC/Fw water were measured at 298 Ka nd 1atm at b = 0.01 mol kg À1 using the Einsteinr elation: D = MSD/6t (where MSD is the mean squared displacemento fa toms and t is time). D = 4.87, 4.83 and 4.77 10 À6 cm 2 s À1 for Nd, Gd and Er,r espectively, with an uncertainty of around 4%.I ti si mportant to note that the self-diffusion coefficient of SPC/Fw water was measured to be 2.35 AE 0.08 10 À6 cm 2 s À1 .T he downward trend in D across the lanthanide series,t hough within statistical uncertainties here, and the order of D values are consistentw ith experimental and other theoretical predictions. [19,44] The enthalpies of solvation were determined as DH solv ¼ DH solution À DH water þ DH corr from simulations at 298 K and 1atm for b(Ln 3+ + ) = 0.01 mol kg À1 . DH water was calculated from a5ns simulation of 5550 SPC/Fww ater molecules at standard temperature and pressure. DH corr = DH B + DH cl + DH comp provides corrections whicha re necessary to compare computed values to experimental data. The formulae for these corrections are provided by Kastenholz and Hünenberger in Ref. [94]. DH B = À319.5 kJ mol À1 and accountsf or errors in solvent polarisation from the use of finite systems with periodic boundaries. DH C1 = À225.8 kJ mol À1 arises due to errors associated with the scheme used to calculate the electrostatic potential. Finally, DH comp = À2.3 kJ mol À1 is the compression work associated with solvation.
DH solv values from C-MD simulations are providedi nT able 2. The experimental values in Ta ble 2h ave been recalculated given that the value used in the originalw ork by Marcus [95] for the hydration enthalpy of ap roton( À1094 kJ mol À1 )h as more recently been reassessed [96] as À1150 kJ mol À1 .F or all cations, the simulations compare reasonably welltothe experiments. [95] DH solv values as predicted by Marcus for Ln 3+ + show an ear linear decrease as af unctiono fa tomicn umber.I no ur work, we find that DH solv ½Gd 3þ is less negative in energy than DH solv ½Nd 3þ ,b ut the generalt rend is for more negative DH solv across the lanthanide series.W ea lso note that the uncertainties in our data are as high as 20 %oft he mean values. The deviationsf rom experiment in the calculated DH solv here for Nd, Gd and Er are around 1, 5a nd 5%,r espectively.T he magnitude of the deviations are not unreasonable when compared to models whichc ontain explicit polarisability.T he force  [44] appearst op erform very well acrosst he lanthanide series with changes < 100 kJ mol À1 when compared with the values of Marcus. [95] Thermodynamic correctionst o DH solv were made but these are not explicitly given by the authors. DH solv was calculated for Nd andG d( see Ta ble2)b yV illa et al. [43] Although their data appear to match well to the experimentalv alues, they applieda+ +211kJmol À1 correction to DH solv to account for ions crossing an interface from the gas phase to the liquid phase. However, Marcus [95] states that the absolute values provided by the experimentalm odelling neglect this effect to properly compares toichiometric quantities to thermodynamic observables. When this correction is removed from DH solv ,t he deviation in the model by Villa et al. from the recalculated experimental values in Ta ble2 is up to 7% (9 %w hen the original hydration enthalpy for ap roton is used).
The C-MD resultsd emonstrate good agreement for our force field with experimental and other theoretical results. Given that our approach to the force-field fitting hasb een systematic, and therefore is easily repeatable for other Ln cations, this procedure provides ar obust andc onsistent tool for exploring thesesystems.

Ion association
We performed two-dimensional US calculations to understand chloridea ssociation to Ln III .H ere, we restrainedb oth the Ln-Cl distance and the number of water molecules in the first Ln coordination shell as defined by S LnÀO in Equation (4). This was necessary given the slow dynamics of water exchange around the trivalent cations.S imulations were performed at b(LnCl 3 ) = 0.01 mol kg À1 in C-MD. Twoa dditional Cl À werer estrained beyondt he second shell of Ln 3+ + to ensure that only single associationb yC l À was sampled. The force constantsfor the additional harmonic restraints were relatively small (k = 10-20 kJ mol À1 ). Te sts showedt hat the additional restraints did not significantly influence the relative changes in free energies between different states in the equilibrium under investigation.
[Nd(OH 2 ) 9 ] 3+ + is thermodynamically stable when Cl À is beyond the second Nd 3+ + coordinations phere, in line with unbiased MD calculations. An energy barrier slightly above thermal energy, k B N A T (k B N A T % 2.478 kJ mol À1 ,w here k B is the Boltzmann constant, N A is theA vogadroconstant and T = 298K), is required forC l À to entert he second Nd 3+ + shella nd form aS ShIP,s ee Figure 4p oint B ,f rom aS SIP, A .T he mean Nd-Cld istance was 5.3 and the TTP coordination geometry was maintained.
To produce the NdCl 2+ + CIP species, D ,ametastable minimum around 5.3 , C ,a nd approximately4k B N A T higheri n energy than the most stable nona-aquoc omplex, is first visited. The structureo ft he coordination complex at C resembles ad istorted SQAP but the value of S LnÀOw is approximately 8.25, given the change in closest water molecule to cation distances. In this water-exchange reaction, [Nd(OH 2 ) 9 ] 3+ + Ð[Nd(OH 2 ) 8 ] 3+ + + H 2 O, the activation energies of forwarda nd reverse reactions were calculatedt ob ea pproximately 5 k B N A T and k B N A T,r espectively.
In the CIP, D ,awater molecule is lost from the Nd 3+ + inner sphere to form am ore stable [NdCl(OH 2 ) 8 ] 2+ + in which the NdÀ Cl distance was 2.91 .T he topology of the energy surfacei ndicatest hat this is ad issociative-activated exchange pathway: one water molecule in the first coordinations phere is lost beforeaC l À enters. Thus, when [Nd(OH 2 ) 8 ] 3+ + is formed at C , there is ac ompetition between aC l À ion and aw ater molecule to move into the first Nd 3+ + shell. The corresponding activation energies were calculated as approx. 2 k B N A T and k B N A T,r espectively.G iven that only very limited water exchange was observed in unbiased MD simulations, this may suggestaslow rate for the forward reaction in Equation (1). The activation energy for this reactioni sa tl east 7 k B N A T. It is possible to lose an additional water molecule to form [NdCl(OH 2 ) 7 ] 2+ + ( E ) through an activated process (~2.5 k B N A T), thought his species is less stable by an energy at least equivalent to thermal energy. DW forG d III (shown in Figure 4) confirms the earlier findings from unbiased MD simulations. Twow ide minima, A and B ,f or [Gd(OH 2 ) 9 ] 3+ + and [Gd(OH 2 ) 8 ] 3+ + are observed when Cl À is beyond the first Gd coordinations phere. However,t he minimum for [Gd(OH 2 ) 9 ] 3+ + is located at S LnÀOw % 8.8 due to an increasei nt he GdÀO water distances. At the furthest sampled distances, the higher coordination state is more stable by < 1kJmol À1 and the energy barrierb etween the two states is around k B N A T. Both coordination states should, therefore, be observedi nu nbiased MD simulations,w ith [Gd(OH 2 ) 9 ] 3+ + being the most probables pecies. As Cl À approaches Gd, the energy difference between the two aquo complexes grows slightly.
An energy barriero f~6 k B N A T separates the GdÀCl SShIP from the more stable CIP at C ,[ GdCl(OH 2 ) 7 ] 2+ + ,i nw hich a water molecule is displaced from the inner sphere on Cl À addition. The CIP has a2 .74 cation-anions eparation distance. Given that [Gd(OH 2 ) 8 ] 3+ + is readily accessible when Cl À is in the second coordination sphere, this marks ac hange in the mechanism for ion association cf. Nd III .R ather than losing water before chloride enterst he first coordinations hell in two steps, in the case of Gd, aw ater molecule is lost when the chloride is added to the first cation coordination spherei naconcerted step through an associative-activated exchange mechanism. This change in mechanism can be ascribed to ad ecrease in cation size on movinga cross the lanthanide series.
For the smallest cation, Er III , DW shows am inimum for [ErCl(OH 2 ) 8 ] 3+ + when Cl À is beyondt he first coordination sphere. This is further confirmationf or ag radualc hange in the most stable water-coordination state for Ln III on moving across the series. Only one minimum is found, A ,f or the ErCl 2+ + SShIP.I n the formation of the CIP, B ,o ne inner-sphere water molecule is replaced by ac hloride ion in ac oncerted step leaving [ErCl(OH 2 ) 7 ] 2+ + .S imilarly to Gd III ,i on association follows associative exchange. The energy barrier to CIP formation was 7 k B N A T and the Er-Cl distance in the CIP was 2.65 .I ti s worth noting that the relative stabilityo f[ ErCl(OH 2 ) 6 ] 2+ + at C is similar to that of the most stable speciesa t B .T hat DW for [LnCl(OH 2 ) 7 ] 2+ + ![LnCl(OH 2 ) 6 ] 2+ + % 5 k B N A T and k B N A T forG da nd Er,r espectively,i sa lso indicative of the increasing cation charge density across the lanthanide series.T his shifting in the stabilities of CIPs to lower coordination states may mean that the mechanism for ion association undergoesf urther changes as the end of the lanthanide series is approached.
An umber of theoretical treatments for ion association in solution have been suggested (see the review by Marcus and Hefter,R ef. [97]). In general,t here is as hort andalong-range contribution to the ion-pairing free energy.T he long-range term is governed by electrostatic forces and can be modelled by the Coulombp otential, taking the solventa sacontinuum with known permittivity.G iven that the permittivity of the mediuma ttenuates the long-range attraction between oppositely charged ions, it is crucial that this is accurately modelled in any simulation. At larger distances still, entropic forces will dominate, and dispersed ions are stabilised as infinite dilution is approached.A ts hort ranges, theories for association differ and this is largely due to the transition from asuitable continuum-solvent model to ad iscrete one and in the model for both the structure and shape of ions.
The two-dimensional PMF surfaces (W(r, S)) in Figure4 were projected onto ao ne-dimensional reactionc oordinate (W(r)), namely the Ln-Cl distance, by taking thermodynamic averages in S LnÀO : An entropy correction to W(r), which accounts for the increasing phase-space volume as af unctiono fc ation-anion radial distance, allows for evaluation of changes in the Helmholtz free energies, where r c is ar eference state. We chose r c to be the maximum value in the sampled Ln-Cl distances. Figure 5s hows the resultingf ree-energy profiles. Calibration of the curves was performed sot hat the mean DA was zero in the region r = 15-16 .I ti si mportant to note that the choice of calibration method can be subjective. The separation distance at which the transition between an ion pair and "free" ions dispersed in solution occurs is not well-defined. Bjerrum showed [98] that the minimum in the probability of finding two oppositely charged ions in solutiono ccurs at ac haracteristic distance, q = z + + z À e 2 / 4 pek B T (where z and e are the ion charge and solution permittivity,r espectively) % 2.1 nm (with ar elative permittivity of 80 for SPC/Fw water), well beyond the limit of two-dimensional US calculations. To verify that the interaction energies between ions were suitablys creened by the solvent, we calculated a free-energy profile (by US simulation with 53 2nsw indows) as af unction of just Nd-Cl separation distance to av alue of r % 2.4 nm. From the free-energy curve in Figure S6 (Supporting Information), we find that DA reaches ap lateau around 1.4 nm. Our calibration procedure is therefore valid, and we considerions separatedbya round1 5 to be dissociated. Figure 5. Relative free energies, DA,for the equilibrium given by Equation (1) (where x = 1) for Nd, Gd and Er.The Coulombic particle-particle interactionenergies, U C ,for + +3e and À1e ions in acontinuum with ap ermittivity equal to that of bulk SPC/Fww ater are also plotted. Curveswere shifted such that all energiesatt he furthest Ln-Cl separation distances were set to zero. Vertical dashed linesshow limits of integration for Nd in the calculation of equilibrium constants(seemain text). Labels indicatethe species associatedwith low energy states in different regionsofL n-Cl distance.
The effect of discrete molecular interactions during ion association is evident in the sequence of barriers that must be overcomeo nm oving from free ions to the CIP.M inima for the SSIP,S ShIP and CIP associates can be seen on decreasing Ln-Cl distance in Figure 5, separated by activation barriersw hich are due to solvation spheres around the cation. Only the barrier for CIP formation is substantial (up to around 15 kJ mol À1 for CIP formation from the SShIP);a ctivation energies in the formation of SSIPs and SShIPs are within 2 k B N A T. Given the small energy barriers and favourable free energies of formation,a relativelyl arge population of SSIPsa nd SShIPs would be expected. There is an oticeable difference in the shape of the activation barrier to forming CIPs for Nd when compared with Gd and Er which is due to the different mechanismsf or ion association. The energy levels for all Ln III ion pairs are very similar. Ta ble 3p rovidest he change in relative free energies from free ions to CIPs.G iven that there is an uncertainty of at least k B N A T in the measurement of free energy changes, it is not possible to identify ac ation as forming more stable ion pairs but the trend shows more favourable binding of Cl À to the lighter lanthanides which has been suggested elsewhere. [28] Thea ctivation energies for CIP formationf rom SShIPsf ollows Nd % Er > Gd, which is consistent with the idea that the frequency of water interchange between the first and second cation coordination spheres reaches am aximum in the middle of the lanthanideseries.
Stability constants for ion pairing can be calculated from changes in free energies between the paired state andd ispersed ions in solution.T he pair correlation-function treatment (see the discussion in the SupportingI nformationa nd Ref. [33] for derivation) has been adopted to calculate stabilityc onstants from the energy profiles in Figure 5. In the case of Ln-Cl ion pairing, the stability constanti sa ne quilibrium constanto f the form, where g ip , g + + and g À are the activity coefficients for ion pairs, cations anda nions, respectively, a is the fractiono fa ssociated ion pairs and c is concentration.G iven that a essentially represents aratio of number densities, we can measurethis through an integration of the free energy profiles in Figure 5. Furthermore-given the nature of the calculation-we can assume that the model represents ion associationa ti nfinite dilution; hence,a ctivity coefficients are one throughout and (1Àa) = 0 within the limits of ion association.
In the above equation, the exponential term is the RDF and r 0 and r m aret he limits over whichi onsa re considered to be associated.H ere, we have used thef eatureso ft he free-energyp rofilest og uide thec hoiceo fi ntegration limits. r 0 wast he minimumi nt he free-energyc urves, whereas r m wast he distance at which DA deviated from zero and this was approximately 15 for all systems. Our approach ensures ac onsistent, empirical comparison between the three systemsunder investigation. Because the adopted formalism is only true for systemsi n which ion activitiese qual concentrations, the resulting equilibrium constants may differ from estimates whichi nclude nonideal factors. We have, therefore, calculated K a values using activity coefficients based on Debye-Hückel theory which can be found in the Supporting Information. However,g iven the semiempirical nature of the theory,w ef ocus our results here on measurementst hat assumen othing aboutt he deviation from ideal behaviour that occurs in solutions at finite concentrations. Thisi sn ot unreasonable, given that stabilityc onstants at infinite dilutiona re often quoted by extrapolating data from experiments performed over af inite concentration range.
The log 10 (K a )v alues are provided in Ta ble 3. We note here that ac oncentration correctioni sm ade to ensure units of dm 3 mol À1 .A lthough, in theory,e quilibrium constants are unitless, the constants calculated using Equation (8) do require a concentrationc orrectioni no rder to be consistent with Equation (7). The log 10 (K a )e valuated from simulations are consistently larger than experimentally determined log 10 (b 1 Cl )v alues which are usually lesst han one. It is important to note that this equilibrium constant accountsf or the formation of associates even with many solvation shells between the cation and anion (g(r)w ill be above one even at 14 ion separation distances).T he decrease in log 10 (K a )a safunction of atomic number is consistent with experimental predictions, [28] demonstratingahigherd egree of ion association for the lighter lanthanides. The gradient in the decrease in log 10 (K a )a safunction of atomicn umber herei sÀ0.01 (R 2 = 0.85). The Supporting Information shows that when activity coefficients deviating from one are included in the calculation, we observet he same trend in the measured log 10 (K a )v alues with the same change as af unction of atomicn umber.H owever,l arger log 10 (K a ) valuesa re obtained using this formalism.
The free-energyp rofiles in Figure5 show that SSIP and SShIP speciesf orm spontaneously in solution with very little energetic penalty.I ndeed, log 10 (K)f or the formation of SShIPs from free ions are above one for all cations studied when the limits of integration are r 1 and r m (see Figure5). The only nonnegligible energy barrier in all of the free-energy profiles is for the equilibrium SShIPÐCIP.I ti s, therefore, informativet oc onsider the equilibrium constant for this step (K CIP )w hich can be calculated according to, Table 3. Free-energy changes, DA,a nd stability constants, K a ,t of orming lanthanide-chloride ion pairs from dispersedi ons in solution and contact ion pairs from solvent-shared ion pairs (K CIP ;see text for details). Statistical uncertainties of around 0.2 kJ mol À1 apply to energy changes, but uncertainties of at least k B N A T = 2.478 kJ mol À1 apply to the data.

DA
log 10 (K a )log 10 where the integral limits refer to ion-separation distances indicated in Figure 5. Ta ble 3g ives log 10 (K CIP )d ata which are within the experimental range for log 10 (b 1 Cl ). [28] Again, relatively higher concentrations of CIPs should be expected for the lighter lanthanides according to these values and the change in log 10 (K CIP )w ith atomic number is linear (R 2 = 0.995) with ag radient of À0.02. These data suggest that whereas ion association is alwaysf avourable,t here will be as ignificant population of weakly associatedi on pairs in solution. log 10 (K CIP )c lose to zero indicates that the population distribution of species will be significantly dependent on solutionc onditions. Given the small driving force for CIP formation and the non-negligible energy barrier to water removal, it is reasonable to assume that significantp opulations of CIPs will only be found at relatively high concentrations of free ions. If we consider the equilibrium between weaklyb ound states including both SShIPs and SSIPs transforming to CIPs, then log 10 (K)v alues are actually negative (À0.646, À0.708 and À0.839 for Nd 3+ + ,G d 3+ + and Er 3+ + , respectively). Although CIPs are energeticallym ore stable than SSIPs and SShIPs, the conformational freedom, and therefore increased entropy,t hat is found for these more weakly bound states makes them collectively more probable in solution. It is likely then, that thermodynamics and kinetics both play ar ole in determining the concentrations of CIPs in solution.

LnCl 3 speciation
Three-dimensionalm etadynamics calculations were performed for LnCl 3 in water using C-MD at 0.1 mol kg À1 .T he collective variabless ampled were S LnÀCl1 , S LnÀOw and S LnÀCl2 ,d efining the coordination between Ln-Cl and Ln-OH 2 in the first Ln III coordination sphere and Ln-Cl within the first two coordination spheres, respectively (see Eq. (4) and Ta ble 2i nt he Supporting Information for detailsa nd parameters). Sampling theseC Vs allows us to compare the relative free energies associated with different[ LnCl x ] (3Àx) species and therefore to rank them in order of thermodynamic stability. In addition, the effect of anions in the second Ln coordination sphere on the stabilities of innersphere complexes can be investigated.
We describe chlorides as being inner sphere,C l i ,w hen occupying the first coordinations phere of Ln III (i.e.,C IPs in the case of ion pairs) ando uter sphere, Cl o ,w hen within the second Ln III coordination sphere (i.e.,S ShIPs in the case of ion pairs) as defined by S LnÀCl1 and S LnÀCl2 ÀS LnÀCl1 ,r espectively.W eh ave chosen to label speciesa ccordingly,f or example, (LnCl i Cl o ) 8 ,w hich refers to a[ LnCl(OH 2 ) 8 ] 2+ + -Cl À speciesw ith one chloride in the Ln III inner sphere along with one solvent-shared chloride and with at otal of eight water molecules immediately surrounding the cation, as indicated by the superscript outside the brackets. Note too, that we have neglected from the labels the total chargeo ft he species.
Plots representing the potential of mean-force energy surface from as ingle NdCl 3 metadynamics calculation are provided in Figure 6. Examples for GdCl 3 andE rCl 3 in water are provided in Figures S7 and S8 (Supporting Information). For clarity, the four-dimensional potential of mean-force surface was projected onto as eries of two-dimensional reaction coordinates. As af unctiono fS LnÀCl1 and S LnÀOw ,w ef ind am inimum for a nona-aquo complex. An energy barrierm ust be crossedb efore chloridea ddition leads to (NdCl i ) 8 ,w hich is the most stable speciesi ns olution, as indicated by the relative energies for some of the sampled species listed in Table 4. It should be noted that whereas Ta ble 4p rovides the statistical uncertainty from averaging multiple metadynamics calculations, there is an uncertainty of at least k B N A T in all of the energy changes reported.
An outer-sphere ion pair,( NdCl o ) 9 ,i ss lightly higheri nenergy than the inner-sphere ion pair and both are more stable than the lanthanide immediatelys urrounded by water.A ll of this is consistentw ith the resultso fo ur US calculations (see Figure 5). The energy differencesb etween these species are smaller than those found in US calculations at b(Ln) = 0.01 mol kg À1 ,h ighlighting the effect of concentration on speciation( also note that entropic corrections for Nd-Cl radial distance are not included here). An inner-sphere complex containingt wo chlorides, ðNdCl i 2 Þ 6 ,w as less stable than Nd immediately surroundedb yn ine water molecules by about k B N A T. The energy barrier to this state from (NdCl i ) 7 was > 10 k B N A T. Gammonse tal. [24] did not detect any [NdCl 2 ] + + in aqueous solu- should predictarelativelys mall population of ðNdCl i 2 Þ 6 based on the calculations presented here. We analysed the structure of ðNdCl i 2 Þ 6 and found ad istorted SQAP geometry with chlorides positioned to maximise the dipolem oment (see Figure S9 in the Supporting Information). In their review,M igdisov et al. [23] discuss the geometries of LnCl 2 + + .I tis suggested that this structure whichw as found for LaCl 2 + + by Petit et al., [99] could resultf rom af avourablec hloride coordination geometry that differs from that of water.T his arrangement does allow for water molecules in the second shell to bind effectivelyt o both water and chloride in the first shell.
An umber of metastable minimaw ere found in the energy landscape which contained two associated chlorides either as ac ombination of inner and outer-sphere, (NdCl i Cl o ) 8 ,o ra s both outer-sphere, ðNdCl o 2 Þ 9 ,i ons. Five minimaa re shown in Figure 6, with ar eactionc oordinate defined by S LnÀCl2 and S LnÀOw ,d ue to ar ange of inner and outer-sphere coordination states and av arying number of water molecules in the first Nd coordination spheres.T hese were all 0 2 k B N A T higheri n energy than the most stable species. (NdCl i Cl o ) 8 was relatively accessible, being separated from (NdCl o ) 8 by an energy barrier of around5k B N A T,i ndicating that the thermodynamic barrier to the formation of aC IP is lowered when Cl À is in the second coordination sphere of the cation. This is an important observation because reduced thermodynamic energy barriers should result in an increasei nt he rate at which cation binding to other ligandst akes place, assuming the mechanismsf or associationa re the same. The energy barriert om ove one chloride from the inner to outer sphere of (NdCl i Cl o ) 8 to form ðNdCl o 2 Þ 9 is around 10 k B N A T. All of the coordination states in which three chlorides were bound to Nd, either in the inner or outerc oordination sphere, were relatively high in energy.E nergyb arriers to adding chloride into coordination complexes which already contained two chlorides weren ot particularly large;h owever,t he disturbance of water structure when three chloridess urround Nd is clearly unfavourable. The relative energy of ðNdCl i 3 Þ 4 compared with the mosts table complex was around + +20 kJ mol À1 .I ti sv ery unlikely,t herefore, that this speciesw ould be observed under standard conditions. The speciation of Gd in chloride solutionsu nder the examined conditions was different to that for Nd. Greaterc hloride coordination to the cation was thermodynamically favourable, with ðGdCl i 2 Þ 5 being the most stable speciesi na ll metadynamics calculations. The complex has ap entagonal-bipyramidal geometry (see the Supporting Information, Figure S9). This structure is consistent with those expectedf rom valence shell electron pair repulsion (VSEPR)t heory.T he energy barriert of orming ðGdCl i 2 Þ 5 from aC IP was, however, above 13 k B N A T,a nd given that there is an on-negligible activation energy to ion pair formation,t he concentration of these species is likely to be extremely low despite their stability. Furthermore, the energy difference betweent his speciesa nd ion pairs is within the range of thermalf luctuations, and so there is not as ignificant driving force for the addition of another anion.A sw as found for Nd, inner-sphere ion pairs, outer-sphere ion pairs and the solvated Gd-cation energies follow the same order predicted in our US calculations, followed by ion-trimer associates with both inner and outer sphere chlorides. Ion quartets with chlorides in the first and second shell werea tl east 5 k B N A T less stable than the most stable species and given that the energy barriers to forming thesea ssociates was above 12 k B N A T,they are unlikely to be detected at 298 K.
The smallest cation (Er) showed the most stable high-coordination states. Ta ble 4s hows that Er bound to two and three chloridesi nt he first coordination sphere were the most stable species. Figure S9, Supporting Information shows that ðErCl i 3 Þ 3 has ad istorted octahedral geometry.A si nt he case of ðGdCl i 2 Þ 5 ,t he geometry of ðErCl i 3 Þ 3 is what can be expected from VSEPR theory.A gain, the activation energies to forming these high-coordination species was considerable:a dding a chloridet ot he inner sphere of ðErCl i 2 Þ 5 required > 40 kJ mol À1 . Following these species, inner and outer-sphere ion pairs and the solvated cation weree quallys table within uncertainties, providing further indication that concentration changes can shift equilibria. Associates containing multiple chlorides in a combination of inner and outer sphere geometries were relatively high in energy,t houghw hen compared to the energies of CIPs, the energy differencesw ere approximately the same as for Gd and Nd. Outer-sphere ion trimers and quartets were unstablebyup to around 30 kJ mol À1 .
The data in Table 4m ay appear to suggestt hat high levels of cation-anion coordination should be expectedi nL nCl 3 aqueous solutions. As stated above, however,t here are considerable activation energies to forming these speciesand so concentrations are likelyt ob ev ery low; nonetheless, the data do show that high-coordination states are either of similarstability to or greater stability than free ions and ion pairs (undert he chosen conditions). Thisi sn ot unreasonablec onsidering the enthalpic gain in couplingp ositive and negative ions, as well as the entropic gain associated with releasing tightly bound water molecules from the cation coordination sphere. It appears that as cation size decreases, the relative stabilities of ion trimers and quartetsi ncreases, and this could be due to the increasingc harge density of the lanthanides acrosst he series.I n general though, the order of stabilities from ion pairs to larger associates is consistent for all three cations. Although it may appear that the range of relative stabilities for different associates increases from Nd to Gd and Er,t he energy differences for for example, adding ac hloride to the outer sphere of the most stable CIP is~2 k B N A T for all cations.C urrent thermodynamicmodelling packages rarely consider species beyond free ions, CIPs and inner-sphere-coordinated ion trimers. For most applications this is probably adequate, particularly because association constantsi nt hese codes can captureadistribution of coordination complex types into as ingle value.H owever,f or a detaileda nalysiso ft he chemistries of the lanthanides in solution, our resultss uggest that consideration of species beyond simple pairs and trimers is important. This is especially the case in hydrothermal fluids in which increasing temperatures should reduce the barrierst oa nion binding.

Conclusions
The force field developed in this work performswell in predicting the behaviour of the lanthanides in water.T his was confirmed by calculations at the DFT level and by comparison to experimentald ata. Furthermore, the enthalpies of solvation highlight that the force field can accurately model the thermodynamics of cationsi nw ater.T he largest deviations in the comparison between calculated and experimental data was for erbium. This is not surprising, given that the force-field fitting procedure relies upon well characterised, stable crystal structures. This computationally inexpensive force field, nonetheless, performs well enough to understand lanthanides in solution. Our methodi sa lso highly systematic, allowing for the addition of other Ln cations and impurities/additives. Two-dimensional US calculations depict clearly that there is at ransitioni nt he relative stabilities of nine versus eight-coordinated Ln III aquo complexes, as has been shown in other studies. [19] Given that the relative free energies for [Gd(OH 2 ) 9 ] 3+ + and [Gd(OH 2 ) 8 ] 3+ + are approximately equal, it is likely that the transition occursa to rc lose to this cation in the lanthanide series. The water-coordination number linearly correlatesw ith the ionic radii of the lanthanides investigated (CN = À0.075Z + 13.5, where Z is atomic number; R 2 = 1). These factorsa lso control the mechanism for ligand exchange. Ac lear shift from dissociated exchange (Nd), akin to an elimination-style process in which the first step is water removal from the inner Ln III sphere, to associatede xchange (Gd and Er), in which water is lost and chloride is added to the inner sphere in ac oncerted step, was apparent from potentialo fm ean-force maps when both Ln-Cl distance and Ln-water coordination were examined. This is consistentw ith ac hange in the mechanism of water interchange, between the first two cation-solvation spheres, aroundt he middle of the series as proposed elsewhere. [15,19] The free-energy profiles in Figure 5s how that, for all of the lanthanides examined, CIPs are the most thermodynamically stable species in the equilibrium in Equation (1) when x = 1. From dispersed ions in solution, the formation of weakly boundS SIPs and SShIPs is both favourable and incurs little energetic cost with thermodynamic barriers within~2 k B N A T. The removal of water to form CIPs from SShIP states incurs an energetic cost of around 15 kJ mol À1 for all cations.T his barrier is very unlikely to be surmounted on the timescales of the equilibrium simulations. Furthermore, CIPs may be inaccessible in some experiments,b ut this energy barrier is certainly accessible on geological timescales andi ne xperiments that are allowed to establish at rue equilibrium.D uvail et al. [42] performed one-dimensional PMF calculations for the binding of Cl to Nd III . They found ab arriert of orming aC IP from aS ShIP of 15 k B N A T which is larger than our estimate of 8 k B N A T. In addition, the SShIP showeda pproximately the same energy as the CIP.H owever,n oe ntropyc orrection was added to the energy profiles to properly consider thermodynamic activation barriers and species stabilities. We believe that if these were included then the Duvail et al. free energies would show the same ordering of thermodynamic stabilities for ion pairs that we present and smaller activation energies to forming ion pairs.
By considering different types of ion pairs on the pathway from free ions to CIPs, our analysis shows that whereas ion association is alwaysf avourable, the equilibrium constantf or the formation of CIPs from SShIPs( K CIP )i sc lose to one for all cations within thermal fluctuations (with at rend to smaller values acrosst he series). The energy change to formingC IPs, calculated using DA = Àk B N A Tln(K), is within the range À0.14 to À1kJmol À1 .T his means that the presence of CIPs will be very dependentu pon the solution conditions. It is interesting to note that our estimate for this equilibrium constant is close to the experimentally determined values but that our measurementso fl og 10 (K a ), which account for the formationo fa ll ion associates,i sm uch larger than the reportedv alues of log 10 (b 1 Cl ), even when activity coefficients deviatingfrom unity are considered (see the Supporting Information). Nonetheless, the trend in decreasing log 10 (K a )a cross the lanthanide series was evident. Our analysisa lso shows that the concentrationso fc ontact pairs should be much lower than the more weakly boundi on pairs (i.e.,S ShIPS, SSIPsa nd beyond) if one can consider a single equilibrium between these types of ion associates: log 10 (K)w as below zerof or all cations and there was at rend to more negative values for the heavier lanthanides.
Previous experimentala nd theoretical studies have failed to providec onsistent conclusions aboutthe nature of chloride association with lanthanides in solution.S tability constantsh ave been calculated [24,25,28] and, whereas these are widely varying, they suggest that ion pairs are stable with respect to dispersed ions. In contrast, spectroscopic measurements and theoretical studies [27,31,32,48,100] suggest that chloride preferentially forms outer-sphere complexes with lanthanides-only at high salt concentrations are contact ion pairs formed. Our study unites these depictions of ion pairing. Ion associationi sa lwaysf avourable, which is not surprising given the strong Coulombic attraction between + +3e and À1e ions in water,b ut weakly bound ion pairs are likely to dominate the equilibrium distribution under relatively mild conditions both thermodynamically and because of the relativelyl arge energy barrierf or the removal of water in the first coordination sphere of Ln III .A t higher concentrations, disturbance in the structure of watersurrounding cationsi sl ikely to lead to ad ecrease in the thermodynamic energy barriers and the equilibrium distribution of ion pairs increases simply because water liberation becomes more favourable. Crucially then, contact-ion-pair formation is controlled both by kinetic andt hermodynamic factors.
Calculationsw hich sampled multiple [LnCl x (OH 2 ) y ] (3Àx) species highlighted the wide number of possible ion associates that can form in solution. Often in speciation analysis, one or two equilibria are considered for the association of chloride to cations;h owever,o ur data show that there is am ultitudeo fe quilibria associated with the formationo fb oth inner and outersphere complexes. Although multiple equilibria can be averaged into just one stability constant, it is important to recognise that estimation of speciesc oncentrations from as ingle stabilityc onstant is not straightforward. Highly coordinated states become more favourable towards the end of the lanthanide series;h owever,t he activation energies to forming these speciesmake them inaccessiblea t298 K. In geological settings, such as hydrothermald eposits, it is likely that there is aw ide range of possible association states that could be considered beyondc ontact ion pairs andt rimers.