Molecular dynamics investigation of water-exchange reactions on lanthanide ions in water/1-ethyl-3-methylimidazolium trifluoromethylsufate ([EMIm][OTf])

Room-temperature ionic liquids are salts that are generally composed of organic cations and inorganic or organic anions. Ionic liquids have attracted attention for a broad range of applications in catalysis,^1–4 organic synthesis,^5,6 separations,⁷ and electrochemistry^8,9 because of their unique properties such as low vapor pressure, non-flammability, and high thermal, chemical, and electrochemical stabilities. These properties can be tuned by choosing an appropriate combination of cations and anions, making ionic liquids task-specific solvents.

Many ionic liquids have been used as solvents for lanthanide solutes because the luminescent,^10–14 magnetic,^15–18 and catalytic^19–22 properties of lanthanides have been extensively exploited for the development of lighting devices, medical diagnosis, catalysis, and organic synthesis. A variety of these properties and applications can be modulated through variation of the lanthanide ions and the coordination environments. Moreover, water-exchange rates influenced by the coordination chemistry of lanthanides are important parameters for characterizing the efficiency of contrast agents for biomedical imaging and diagnosis and catalysts in organic reactions.^23,24 Recently, ionic liquids have become important for the coordination chemistry of lanthanides because lanthanide ions in ionic liquids often exhibit unique properties compared to conventional organic solvents.^14,16,22 A number of properties and behaviors of lanthanides in ionic liquids have been studied experimentally^14,16,22 and theoretically.^25,26 However, to our knowledge, water-exchange kinetics of lanthanide ions in ionic liquids are scarcely discussed, and more comprehensive structural and dynamics studies of these processes in ionic liquids will be helpful for the development of lanthanide-based catalysts for diverse applications in both the chemical and biomedical fields.

Molecular dynamics (MD) simulations can explicitly represent solvents and, therefore, have become powerful tools to study solvation characteristics and the water-exchange processes.^25–33 However, rates of water exchange are influenced by metal–water and water–solvent interactions, so an accurate force field is required for appropriately modeling these intermolecular interactions. The atomic-multipole-optimized-energetics-for-biomolecular-application (AMOEBA)³⁴ force field has been used to simulate water-exchange dynamics for various metal ions because it can accurately describe the interactions between metal ions and water molecules. Recently, our group has developed a protocol to incorporate Gaussian electrostatic model-distributed multipoles (GEM-DMs)³⁵ into the AMOEBA force field and was able to successfully estimate the bulk properties for several ionic liquids.³⁶ In addition, our previous MD study based on the AMOEBA force field was able to reproduce the experimental water-exchange behaviors that revealed that the trend of water-exchange rates in water/[EMIm][EtSO₄] for different lanthanide ions is opposite to that observed in water.^33,37 Our simulation results demonstrated that the dissociation of a water from the first shell is a key step for water-exchange events in water/[EMIm][EtSO₄]. Smaller lanthanide ions undergo faster water exchange in water/[EMIm][EtSO₄] because they bind to the first shell [EtSO₄]⁻ anions more tightly than larger lanthanide ions, resulting in more steric effects on the first shell water molecules and facilitating the dissociation of water.

Lanthanide trifluoromethylsufates (or triflates) have been extensively used as Lewis acids in a wide range of organic-synthesis and catalysis applications in aqueous solution because these acids are tolerant of water.^38–40 In addition, because of the high thermal and electrochemical stability of triflate ionic liquids, they are of interest as the heat transfer materials and the electrolytes in the electrochemical cells.^41,42 However, water-exchange reactions of lanthanide ions in triflate ionic liquids have not been studied, and more comprehensive studies will be important for further applications of these ionic liquids. In this paper, we explore the water exchange processes for different lanthanide ions in water/[EMIm][OTf]. We follow our previous methodology³⁶ to develop the AMOEBA force field parameters for [OTf]⁻, which are not available in the current AMOEBA force field. With the optimized force fields, we can simulate rates of water-exchange reactions in water/[EMIm][OTf]. The solvation characteristics, transport properties, and solvated structures from the MD trajectories were further examined to understand the water-exchange behaviors of lanthanide ions in water/[EMIm][OTf].

The force-field parameters for lanthanide ions and the [EMIm]⁺ cation have been developed previously,³³ and these optimized parameters were used for simulation of water-exchange reactions on lanthanide ions in water/[EMIm][OTf]. Here, we employed the same computational procedures to develop the required AMOEBA force-field parameters for [OTf]⁻.³⁶ Briefly, we optimized the structure for a [OTf]⁻–water dimer at the MP2 level^43–47 with a 6-311G(d,p) basis set using the Gaussian software package.⁴⁸ The relaxed one-electron density from the optimized structure was employed to calculate the distributed multipoles via the Gaussian electrostatic model-distributed multipole (GEM-DM) method.³⁵ At the optimized geometry, second-order symmetry-adapted perturbation theory (SAPT2)^49–51 energy decomposition analysis was performed using the Psi4 package⁵² to determine individual components of intermolecular interactions for different [OTf]⁻–water dimers as a reference. All parameters were optimized by comparison to different quantum-mechanical calculations including optimized structures, intermolecular interactions, and energy decomposition analyses as previously described.³⁵ MD simulations were performed with 216 ion pairs of [EMIm][OTf] to test the accuracy of the developed force field for [EMIm][OTf]. The simulated densities and heats of vaporization were compared with reported experimental values.^53,54 The parameters for the [EMIm][OTf] ion pair are provided in the supplementary material.

Simulations for three lanthanide ions—Gd³⁺, Dy³⁺, and Ho³⁺—were performed with a mixture of water and [EMIm][OTf]. To match the experimental NMR conditions (the volume ratio 1:19 of water to [EMIm][OTf]), 99 water molecules were placed randomly in a box of 180 [EMIm][OTf] ion pairs. Three [EMIm]⁺ cations were replaced with a lanthanide ion for each simulation studied to keep the systems electroneutral. All simulations were performed with the AMBER14 package⁵⁵ and the AMOEBA polarizable force field with the above-determined parameters. The simulations were carried out in an isothermal isobaric ensemble (NPT) at 298 K and 1 bar with a 1 fs time step using an 8.5 Å cutoff for non-bonded interactions and an 8.0 Å cutoff for the smooth particle mesh Ewald (PME) method⁵⁶ for long-range electrostatics. The Berendsen thermostat and barostat⁵⁷ were used to control temperature and pressure with relaxation times of 1 and 3 ps, respectively. All production trajectories were generated for at least 20 ns, and snapshots were recorded every 100 fs. All properties—such as radial distribution function, molecular diffusion coefficient, and mean residence time—were calculated from three independent MD trajectories with the standard deviation used as the uncertainty.

The experimental method and procedures used to measure water-coordination numbers and water-exchange rates of lanthanides in water/ionic liquids have been reported,³⁷ and the data collection in this work followed those procedures with the change of [EtSO₄]⁻ from that report being replaced by [OTf]⁻ in this study. The luminescence-decay measurements of Eu(OTf)₃ and Tb(OTf)₃ (Table S9 of the supplementary material) were used to determine an average water-coordination number of 1.0 for Eu³⁺ or Tb³⁺ in water/[EMIm][OTf] (1:19, v/v). Data from variable-temperature ¹⁷O-NMR measurements of lanthanide triflates [Gd(OTf)₃, Dy(OTf)₃, and Ho(OTf)₃] and their diamagnetic reference [Y(OTf)₃] collected in the experimental determination of water-exchange rates are provided in the supplementary material (Tables S9–S14 and Figs. S10–S12).

Water-exchange rates of three lanthanide ions in water/[EMIm][OTf] were determined using ¹⁷O-NMR relaxation-rate measurements. The results show that water-exchange rates of lanthanide ions in water/[EMIm][OTf] increase as a function of charge density (the ratio of charge to ionic radius). The increase of water exchange rates from Gd³⁺ to Ho³⁺ was also observed in water/[EMIm][EtSO₄]. However, the exchange rate increases more rapidly with the increasing size of metal ions in water/[EMIm][OTf] than in water/[EMIm][EtSO₄]. To understand the water-exchange reactions in this solvent, we started with the development of force-field parameters for [EMIm][OTf]. To validate the force field for [EMIm][OTf], the density, the heat of vaporization, and the self-diffusion coefficient for [EMIm][OTf] were simulated. With the optimized force-field parameters, we simulated the water-exchange dynamics of lanthanides in water/[EMIm][OTf].

Table I lists the experimental and simulated densities and heat of vaporization for [EMIm][OTf] for four different temperatures. A density of 1.413 g/cm³ predicted from simulation at 298 K is 2.1% higher than the experimental values of 1.384 g/cm³ using the original van der Waals (vdW) parameter from the AMOEBA potential. In addition, the calculated heat of vaporization (97.5 kJ/mol) is underestimated by 29.3% compared to the experimental ΔH^vap of 137.9 kJ/mol. Therefore, the van der Waals parameters for [OTf]⁻ were optimized by comparing the calculated density and ΔH^vap with their experimental counterparts at a single temperature (298 K) to improve agreement with the experimental data. After the van der Waals parameters for [OTf]⁻ were optimized, the calculated density (1.394 g/cm³) and heat of vaporization (136.6 kJ/mol) at 298 K show errors below 1% with respect to the experimental counterparts. The calculated density slightly decreases with increasing temperature (Table I) and is in good agreement with experimentally reported data, with errors below 2% for all temperatures.

Radial distribution functions, g(r), were calculated for each lanthanide ion (Gd³⁺, Dy³⁺, and Ho³⁺) with various solvent molecules (water, [OTf]⁻, and [EMIm]⁺) to examine the coordination environment of lanthanide ions in water/[EMIm][OTf]. The peak positions of g(r) indicate the location of solvent coordination shells around the lanthanide ions, and the corresponding integral values of g(r) peaks, CN(r), represent the average coordination numbers of solvent molecules. The peak positions and the corresponding integrated values are summarized in Table II.

The g(r) for the lanthanide ion with inner-sphere water are displayed in Fig. 1 and are given by water O (Ln–O_water g(r), solid line) and H atoms (Ln–H_water g(r), dashed line). In Fig. 1, Ln–O_water g(r) display two well-resolved peaks, indicating the population of water in the first and second coordination spheres of the lanthanide ion. The first sharp peak appears at 2.23 ± 0.01 Å, indicating a population of water in the first coordination sphere of the lanthanide ion. The peak position of Ln–O_water g(r) for all lanthanide ions is located at almost the same distance and shows no dependence on the size of the lanthanide ions. However, the maxima of these peaks are centered at shorter distances compared to the previously reported distances of 2.35–2.38 Å in water and 2.40–2.43 Å in water/[EMIm][EtSO₄].³³ This difference indicates that the binding strength between the lanthanide ions and water in water/[EMIm][OTf] is stronger compared to only water or water/[EMIm][EtSO₄], possibly resulting from less steric repulsion between water and [OTf]⁻ than water–[EtSO₄]⁻ and water–water in the first coordination shell.

The first and second peaks of Ln–O_water g(r) are well separated by a minimum, and the minimum reaches zero. In addition, integration of the first Ln–O_water g(r) gives an integer of 1.0. These results indicate that only one water molecule stays in the first shell and no water-exchange processes occur between the first and second shells during the simulations.

The second Ln–O_water g(r) peak is broader than the first one and appears at 4.47 Å for Gd³⁺, 4.51 Å for Dy³⁺, and 4.55 Å for Ho³⁺. The second and third solvation shells are also separated by a minimum, but the minimum after the second peak does not reach the x-axis, indicating that second-shell water exchanges with outer-sphere water. The integral values of the second Ln–O_water g(r) peaks show that the number of water molecules in the second coordination shell is between one and two due to water-exchange events between the second shell and the bulk.

Similar to Ln–O_water g(r), Ln–H_water g(r) display two well-defined peaks corresponding to the first and second solvation shells. The first peak of Ln–H_water g(r) has zero overlap with the first one of Ln–O_water g(r), indicating that the water dipole in the first coordination shell is strongly oriented toward the lanthanide ion. The second Ln–H_water g(r) peak appears at 5.0–5.1 Å and is partially overlapped with the second Ln–O_water g(r) peak, indicating that the second shell water molecules are also strongly oriented by the lanthanide ion but have more flexibility compared to the first shell water molecules. The higher flexibility of the second shell water is consistent with water-exchange with the bulk.

The Ln–O_OTf g(r) for the lanthanide ions and the O atoms of [OTf]⁻ anions for the distribution of [OTf]⁻ anions were also calculated (Fig. 2, left panel). For all lanthanide ions, the Ln–O_OTf g(r) patterns are similar, and a well-resolved peak at the distance of 2.7–2.8 Å indicates the location of the first coordination shell formed by the O atoms of [OTf]⁻ anions. The first peak of Ln³⁺–O_OTf g(r) appears at 2.79 Å for Gd³⁺, 2.77 Å for Dy³⁺, and 2.75 Å for Ho³⁺. Integration of the first Ln–O_OTf g(r) peaks show that the lanthanide ions are coordinated with nine O atoms from six [OTf]⁻ anions in the first shell. In addition, the first peak of Ln–O_OTf g(r) is broad compared to that of Ln–O_water g(r), and the broader peak for Ln–O_OTf g(r) is partly due to interconversion of coordination modes of [OTf]⁻ between bidentate and monodentate during the simulation.

The Ln–[EMIm]⁺ g(r) show zero population of [EMIm]⁺ cations in the first coordination shell because of the repulsion between the lanthanide and [EMIm]⁺ cations (Fig. 2, right panel). When the distance of the lanthanide ions and [EMIm]⁺ is longer than 4 Å, the Ln–[EMIm]⁺ g(r) become non-zero, indicating that [EMIm]⁺ cations occupy the second coordination shell of lanthanide ions.

Determination of diffusion coefficients can provide useful information about the molecular mobility and kinetics of water-exchange reactions in a given solvent environment. The molecular diffusion coefficients (D) can be determined by using Einstein’s relation⁶¹ 

where MSD is the mean-square displacement of molecules over a period of time (t). In the diffusive regime, the MSD grows linearly with time, which enables calculation of diffusion coefficients from the slopes of the plots of MSD as a function of time. The molecular diffusion coefficients (D_H2O, D_OTf, and D_Ln, Ln = Gd³⁺, Dy³⁺, or Ho³⁺) in water/[EMIm][OTf] were calculated from the slopes of MSD plots from 500 to 15 000 ps (Table III).

To study the effect of the metal ion on molecular diffusion, the diffusion coefficients of water and lanthanide ions were determined for Gd³⁺, Dy³⁺, and Ho³⁺ in water/[EMIm][OTf] (Table III). All three lanthanide ions have similar diffusion coefficients in water/[EMIm][OTf]: 1.53 × 10⁻¹⁰ cm²/s for Gd³⁺, 2.88 × 10⁻¹⁰ cm²/s for Dy³⁺, and 2.35 × 10⁻¹⁰ cm²/s for Ho³⁺. The effect of the metal ion on the diffusion coefficients of water molecules in water/[EMIm][OTf] is also not significant, with the diffusion coefficients of water (D_H2O) corresponding to 1.94 × 10⁻⁸ in Gd³⁺/water/[EMIm][OTf], 1.80 × 10⁻⁸ cm²/s in Dy³⁺/water/[EMIm][OTf], and 1.33 × 10⁻⁸ cm²/s in Ho³⁺/water/[EMIm][OTf].

The calculated diffusion coefficients of water (D_H2O) were compared to our previously reported values in water or water/[EMIm][EtSO₄] to investigate the solvent effect on diffusion of water molecules. Table III shows that the diffusion coefficients of water in water/[EMIm][OTf] are in the range of 1.33 × 10⁻⁸ to 1.94 × 10⁻⁸ cm²/s, slower than those in water (1.59 × 10⁻⁵ to 1.63 × 10⁻⁵ cm²/s).³³ The slower molecular diffusion of water in water/[EMIm][OTf] is primarily due to the formation of intermolecular hydrogen bonds between water and [OTf]⁻ anions. The diffusion coefficients of [OTf]⁻ are in the range of 1.45 × 10⁻⁹ to 2.39 × 10⁻⁹ cm²/s in water/[EMIm][OTf], suggesting that the hydrogen bonds between water and [OTf]⁻ anions slow the diffusion of water in water/[EMIm][OTf]. The diffusion of molecules influenced by hydrogen-bond interactions was also observed in our previous study in water/[EMIm][EtSO₄].³³ In addition, water molecules diffuse more rapidly in water/[EMIm][OTf] (1.33 × 10⁻⁸ to 1.94 × 10⁻⁸ cm²/s) than in water/[EMIm][EtSO₄] (1.40 × 10⁻⁸ to 8.64 × 10⁻⁹ cm²/s), and this difference is consistent with the ¹⁷O-NMR experimental results demonstrating that water-exchange rates of lanthanide ions are faster in water/[EMIm][OTf] than in water/[EMIm][EtSO₄].

Two methods can be used to calculate the water-exchange rates, namely, the survival function proposed by Impey et al.⁶² and the direct method,⁶³ accounting for all incoming and outgoing water molecules throughout the whole simulation trajectory. In this work, the direct method was chosen because this method treats the associative and dissociative exchange processes equally. Similar to the Impey procedure, the direct method also needs the time parameter t* for defining a real exchange event. In this work, the parameter t* is taken to be 2 ps, which means that if a water molecule enters or leaves the second shell (i.e., 4.0 Å < distance of Ln–O_water < 5.5 Å) for a time period longer than 2 ps, the event is counted as one real exchange process. The average residence time (τ_ex) of water can be defined as

where t_sim is the whole simulation time, CN_av is the average number of water molecules in the second shell, and N_ex is the number of water-exchange events between the second shell and the bulk. The distances of the lanthanide ion with water O atoms at various solvation shells were analyzed throughout the whole simulation time to record all exchange events of water. In Fig. 3, the Ln–O_water distance trajectories show that no water exchange events occurred from the first solvation shell of lanthanide ions. Nevertheless, a number of second shell water-exchange processes were observed during the MD simulations. The calculated water-exchange rates (k_ex), the inverse of mean residence time (τ_ex), were determined from the simulation time of 20 ns and listed in Table V. As shown in Table V and Fig. 3, the calculated water-exchange rates are in the order of Gd³⁺ < Dy³⁺ < Ho³⁺. The increasing trend for the rates is similar to that in water/[EMIm][EtSO₄] and is consistent with the ¹⁷O-NMR experiments (Tables IV and V).

Further information regarding water-exchange events can be obtained by examining the snapshots of lanthanide ions with the inner-sphere solvent molecules from the MD trajectories. Figure 4 shows a representative example of the water-exchange process of a lanthanide ion in water/[EMIm][OTf]. At the beginning of the simulations (t = 0 ps in Fig. 4), all lanthanide ions were coordinated by six [OTf]⁻ anions and one water molecule in the first solvation shell. Six [OTf]⁻ anions form a trigonal prism structure, and a water molecule is located at one of the square faces of the trigonal prism structure. The second solvation shell consists of two water molecules and two [OTf]⁻ anions. Two water molecules form hydrogen bonds with two H atoms of the first-shell water, and the O atoms of the second-shell water molecules are pointed toward the lanthanide ion, showing strong dipole interactions between the water O atoms and the lanthanide ion. During the simulation, both the [OTf]⁻ anions and water molecules undergo rapid rotational motions around the lanthanide ion. In some occasions, the rotational motions cause the second shell water to form hydrogen bonds with the O or F atoms of the first shell [OTf]⁻, disrupting the interaction of the lanthanide ion with the second shell water molecules and leading to the dissociation of a water from the second shell to the bulk (blue water at t = 7706–7707 ps in Fig. 4). Once the water leaves the second shell, the outgoing water becomes randomly oriented toward the lanthanide ion. The entering water can be the original outgoing water or other water molecules from the outer-sphere.

Intermolecular interaction energies of Ln–[OTf]⁻ and Ln–water were calculated for selected snapshots from the simulation trajectories to understand the effect of the metal ion on the water exchange rates in water/[EMIm][OTf]. All intermolecular interactions were determined using the AMOEBA force field and TINKER program.⁶⁴ The averaged interaction energies between the lanthanide ion and [OTf]⁻ anions and between the lanthanide ion and water molecules are summarized in Fig. 5. The interaction energies of the lanthanide ion with [OTf]⁻ anions and water are dependent on the size of the lanthanide ion. In the first solvation shell [Fig. 5(a)], both [OTf]⁻ anion and water bind more tightly to Ho³⁺ than Gd³⁺ due to the larger charge density of Ho³⁺. On the other hand, interaction energies of the lanthanide ions with the second shell [OTf]⁻ anions and water become less negative when the size of the lanthanide ion decreases from Gd³⁺ to Ho³⁺[Fig. 5(b)]. A possible cause for the decreasing interactions of Ln–water and Ln–[OTf]⁻ involves a larger screening effect resulting from the stronger interaction between the smaller lanthanide ion and the first shell [OTf]⁻. The first shell [OTf]⁻ anions bind to Ho³⁺ more strongly than Gd³⁺, making it more difficult for Ho³⁺ to interact with the second shell [OTf]⁻ and water compared to Gd³⁺. As a result, the second shell water molecules around Ho³⁺ are more labile than those around Gd³⁺, and the trend of Gd³⁺ < Dy³⁺ < Ho³⁺ is seen for water-exchange rates in water/[EMIm][OTf].

Experimental results show that the water-exchange rates of lanthanide ions in water/[EMIm][OTf] increase in the order of Gd³⁺ < Dy³⁺ < Ho³⁺. The trend is similar to that observed in water/[EMIm][EtSO₄]. To gain further insight into the water-exchange events in water/[EMIm][OTf], we performed classical MD simulations using the AMOEBA polarizable force field. Our simulations are in good agreement with experiment with respect to the trend of water-exchange rates. By analyzing the MD trajectories, we observed that the water-exchange event for a lanthanide ion occurred between the second shell and bulk in water/[EMIm][OTf]. The exchange rate of the second-shell water with the bulk water depends on the relative binding strength of the lanthanide ion with the first shell [OTf]⁻ anions. The [OTf]⁻ anion binds to the smaller lanthanide ions more tightly, but this causes the more screening effect to the second shell water molecules. The stronger screening effect weakens the interaction of the smaller lanthanide ion with the second shell water, facilitating water exchange from the second shell of the smaller lanthanide ion. As a result, the trend of Gd³⁺ < Dy³⁺ < Ho³⁺ is seen for water-exchange rates in water/[EMIm][OTf]. Although the trend of water-exchange rates in water/[EMIm][OTf] is similar to that in water/[EMIm][EtSO₄], our simulation indicates that the process of water exchange for lanthanide ions in water/[EMIm][OTf] is different from that in water/[EMIm][EtSO₄]. This study suggests that the second-shell water exchange should be considered when interpreting ¹⁷O-NMR data, and it provides a possible mechanism for manipulation of water-exchange dynamics through selection of different ionic liquids as solvents.

See supplementary material for definitions for all local frames, multipoles, and intra- and inter-molecular interactions for the [EMIm][OTf] ion pair. MD simulation results for density and heat of vaporization, self-diffusion coefficients, and mean residence time.

Y.-J. Tu thanks Wayne State University for a Thomas Rumble Fellowship. M.J.A. acknowledges the National Science Foundation for support (Grant No. CHE-0955000). The authors also thank Wayne State computing grid for computational resources and R. L. Chamika Lenora for helpful discussion.