Current Status of AMOEBA–IL: A Multipolar/Polarizable Force Field for Ionic Liquids

Computational simulations of ionic liquid solutions have become a useful tool to investigate various physical, chemical and catalytic properties of systems involving these solvents. Classical molecular dynamics and hybrid quantum mechanical/molecular mechanical (QM/MM) calculations of IL systems have provided significant insights at the atomic level. Here, we present a review of the development and application of the multipolar and polarizable force field AMOEBA for ionic liquid systems, termed AMOEBA–IL. The parametrization approach for AMOEBA–IL relies on the reproduction of total quantum mechanical (QM) intermolecular interaction energies and QM energy decomposition analysis. This approach has been used to develop parameters for imidazolium– and pyrrolidinium–based ILs coupled with various inorganic anions. AMOEBA–IL has been used to investigate and predict the properties of a variety of systems including neat ILs and IL mixtures, water exchange reactions on lanthanide ions in IL mixtures, IL–based liquid–liquid extraction, and effects of ILs on an aniline protection reaction.


Introduction
The study of ionic liquid (IL) solutions by means of molecular dynamics (MD) or Monte Carlo (MC) computational simulations have become a useful tool to study these systems. These approaches can provide significant insights on structural, thermodynamic and transport properties. Many of these approaches employ classical force fields (FFs), which employ bonded (bond length, angle, etc.) and non-bonded (Coulomb, Van der Waals) terms to approximate the energy of these systems [1].
Several groups have developed pairwise additive force fields (FFs) for the simulation of ILs [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. A large number of studies have been devoted to the modeling of ILs by MD based on these FFs [10,. In general, trends for thermodynamic, structural and transport properties are reproduced. However, the limited accuracy of current FFs results in over-or under-estimation of some of these properties. Case in point, errors in calculated enthalpies of vaporization compared to experiment may be as high as 50% [11,37]. One more issue that needs to be taken into account is the high viscosity of these liquids that results in the need for long simulation times to obtain converged data for transport properties.
The parametrization philosophy of AMOEBA-IL is based on the accurate reproduction of quantum mechanical (QM) data for dimer and oligomer systems, as well as the reproduction of bulk properties. In all cases, the bonded terms have been taken from the original AMOEBA parameters where available. For the bonded AMOEBA parameters that have not been previously reported, torsional scans have been obtained at the MP2/6-311G(d,p) level [102].
All parameters for the non-bonded terms for AMOEBA-IL are fitted using counterpoise corrected QM inter-molecular interaction energies and QM energy decomposition analysis (EDA) for representative dimers [102][103][104][105][106]. The QM EDA data is employed to fit and compare the individual Coulomb, polarization and Van der Waals AMOEBA terms to reduce the error of each term. The original parametrization for the dimethyl imidazolium-based systems employed the restricted variational space (RVS) approach for the QM EDA reference using dimers in a single orientation at different distances [102,103].
For the most recent parametrizations, the symmetry adapted perturbation theory (SAPT) method has been employed to obtain the reference components, coupled with randomly oriented molecules at different distances (see Figure 1) [105,106]. This newer approach takes advantage of an efficient implementation of SAPT in Psi4 [107], which enables the use of higher accuracy QM reference, as well as better representation of the dimer surface using randomly oriented dimers.
As a specific example, one of the cations included in AMOEBA-IL is spyrocyclic pyrrolidinum [SPyr] ( Figure 2). This cation comprises two five-membered rings fused by the quaternary ammonium atom. The two cycles are compsed of sp 3 C atoms, which allow bending of the two rings with respect to the central N. The parametrization of the multipoles and VdW for this cation was performed using a set of thirty randomly-oriented [SPyr]-water dimers and validated against seventy [SPyr][BF 4 ] randomly oriented dimers ( Figure 1). Given the [SPyr] internal felxibility, two sets of parameters were developed comprising a set with no internal polarizable group (1G) and a set with three polarizable groups (3G) (Figure 2) to investigate the role of intra-molecular polarization (see below). The comparison of each non-bonded term with respect to QM EDA data provides a way to determine how each individual contribution is reproduced by AMOEBA-IL. This comparison allows a better understanding of the performance of the individual components of the force field.  Distributed atomic multipoles for AMOEBA can be obtained using the Gaussian distributed multipole analysis (GDMA) approach from Stone [108]. Alternatively, we have shown that strictly convergent distributed multipoles can be obtained from Hermite Gaussian functions using the Gaussian electrostatic model (GEM) fitting approach [109][110][111]. The Coulomb and polarization interactions are compared against the corresponding EDA counterparts. Additionally, it is possible to approximate intra-molecular polarization via polarizable groups as introduced by Ren and Ponder and employed for AMOEBA-IL (see Figure 1) [98,105].
The Van der Waals term is fitted by using the energies obtained by subtracting the AMOEBA calculated Coulomb and polarization from the total counterpoise-corrected QM energy [102][103][104][105]. Thus, the Van der Waals term effectively includes not only exchange and dispersion interactions, but also charge transfer effects and other errors in the non-bonded terms such as Coulomb penetration. In some cases it is advantageous to employ bulk properties to further refine the Van der Waals parameters [102,105,112].
This approach has resulted in parameters that provide accurate description of neat ILs and IL solutions (vide infra). For example, for [EMIM][EtSO 4 ] the AMOEBA-IL calculated density and heat of vaporization at 298 K results in errors less than 1% [103], compared with errors ≈1.5% and >5% for density and ∆H vap respectively with conventional point charges [113]. Similarly, for [EMIM][OTf], AMOEBA-IL results for these properties also show errors below 1% [104] compared with ≈1.5% for density and >1% compared with >7% for the heat of vaporization at 298 K [113,114]

Water Exchange Dynamics on Lanthanide Cations
Lanthanide ions (Ln-ions) are employed as contrast agents for biomedical imaging because of their luminescent and magnetic properties. One feature that can affect the contrast agent efficiency is the rate of ligand exchange on the complex, in particular that of water, [115,116], which can be modulated by altering the coordination environment. Thus, understanding the water exchange mechanism in Ln-ions can provide important insight for contrast agent development.
Some of us performed MD simulations employing AMOEBA to study the structure and dynamics on these solvent-exchange processes. The AMOEBA parameters for Gd 3+ , Dy 3+ , and Ho 3+ ions were obtained by comparing energies calculated with the respective AMOEBA parameters with the interaction energies of Ln-water dimers obtained at the MP2/SDD/6-311G(d,p) level of theory (Stuttgart's small core quasi-relativistic effective core potential (SDD) for Gd 3+ , Dy 3+ , and Ho 3+ ions and 6-311G(d,p) basis sets for H and O atoms) and previously reported QM EDA data for different Ln-water dimers [117]. Polarizabilities of Lanthanide trivalent cations have been previously reported by Marjolin et al. [117][118][119].
The parameters for the lanthanide cations were tested by determining the coordination numbers and radial distribution functions (g(r)) in water and in the mixtures with the two ILs, and comparing to experimental results (see Figure 3). Consistent with experimental data, the Ln-Ow distances are seen to decrease as the ionic radius decreases from Gd 3+ to Ho 3+ . Additionally, the AMOEBA-based simulations predict between eight and nine water molecules in the first hydration shell due to water-exchange events between the first and the second shells in agreement with experimental results [103].
In the water/[EMIm][EtSO 4 ] system, Ln-ions can be surrounded by water and [EtSO 4 ] − anions ( Figure 3). The two peaks observed in Figure 3 for the first coordination shell suggest the Ln-ions are all coordinated by four [EtSO 4 ] − anions. Three anions coordinate the respective lanthanide cation in a bidentate fashion, and one more as a monodentate ligand. Additionally, the first coordination shell of the Ln-ions contains one or two water molecules due to water-exchange events, resulting in the two observed peaks.
Water-exchange rates were calculated using the survival function method [125] (for water/[EMIm][EtSO 4 ] system) and the direct method [126] (for water/[EMIm][OTf] system); both methods need a time parameter (t * ) for defining a real exchange event. A water-exchange event relates the time difference between a water molecule coming/leaving from the first solvation shell. In water, water-exchange rates exhibit the trend of Gd 3+ > Dy 3+ > Ho 3+ , but slower in water/[EMIm][EtSO 4 ] and water/[EMIm] [OTf]. The calculated water-exchange rates show the same trends as the 17 O-NMR experiments, albeit they are slightly faster (Table 1). Notably, the water exchange rate (both experimental and calculated) in neat water is observed to decrease as the Ln-ion atomic number increases. Conversely, in both water/IL mixtures this trend is reversed, that is, the water exchange rate decreases with increasing atomic number of the Ln-ion. Water-exchange events between the first coordination shell and the bulk along the MD trajectories were analyzed by measuring the distance between the Ln-ion and the oxygen of the water molecules. (Figures 4-6) in different solvent systems. In water, at the beginning of the simulation, the first hydration shell of Ln 3+ has eight water molecules, forming a square antiprism (SAP) geometry. Along the simulation, one water molecule from the bulk joins the first hydration shell, and a nine-coordinate Ln-aquo complex is formed with a tricapped trigonal prism (TTP) geometry, followed by a water release to the outer hydration shells, and rearranging back to the SAP geometry. Based on their ionic radii (Gd 3+ > Dy 3+ > Ho 3+ ), smaller Ln-ions form tighter aquo complexes with the neighboring water molecules, affecting the accessibility for the incoming water and thus impacting to the water-exchange rates. Therefore, our simulations suggest that the water exchange in neat water follows an associative mechanism based on the overlap of Ln-water distance trajectories for the incoming and outgoing water molecules in the first hydration shell.
For the water/[EMIm][EtSO 4 ] system, initially the Ln 3+ ion is coordinated with nine O atoms from four [EtSO 4 ] − and two water molecules. No exchanges of water or other ligands occurred in the first nanoseconds of MD simulations. Interestingly, [EtSO 4 ] − experienced rapid spin motions around the Ln-ion, resulting in the occasional increase of the Ln-water distance (dissociation of a water molecule), losing a water molecule and promoting a water-exchange process. The residence time trend of a water in the first shell is opposite to the one for water-exchange rates in water (Ho 3+ < Dy 3+ < Gd 3+ ) and depends on the charge density of the Ln-ion. [EtSO 4 ] − anions bind strongly to smaller Ln-ions, increasing the steric effects in the first coordination shell, impeding the water/Ln 3+ binding and resulting in a more facile release of water molecules from Ln 3+ . No water-exchange event was observed in all MD trajectories when two first shell water molecules were adjacent to each other at the beginning of the simulation. Based on these results, the water exchange mechanism for the water/[EMIm][EtSO 4 ] mixture corresponds to a dissociative mechanism in contrast to the neat water mechanism.
For the water-exchange process in water/

Liquid-liquid Extraction of Benzene from Dodecane-Benzene Mixture
Ionic liquids are attractive solvents for liquid-liquid extraction due to their unique properties (low vapor pressure, reusability, thermal and chemical stability). Several researchers have obtained experimental evidence of the extraction of benzene (PhH) and other aromatic compounds from hydrocarbon mixtures like gasoline using ILs [127][128][129][130][131][132][133]. We have recently employed AMOEBA-IL to investigate the extraction of PhH from a gasoline-model using two imidazolium-based ILs [106].
Two systems consisting of a gasoline-model  (Figure 8). The density profiles along the z direction were used to measure the PhH extracting capabilities, and the spatial distribution functions (SDF) were employed to gain further insights on the interactions within the studied mixtures.   (Figure 10a,c), show that at short range, there are stacking interactions between aromatic groups, and these interactions are similar in binary and ternary systems (Figure 10b,d). On the other hand, at long range, there are more PhH/[DMIM] + interactions compared with PhH/[EMIM] + . These differences arise from the imidazolium ion's alkyl chain length. 4 ] ternary system both along the edge and plane of the ring (Figure 11b,d).

Computational/Experimental Characterization of Spyrocyclic Pyrrolidinium/Tetrafluoroborate [Spyr][BF −
4 ] ILs have been studied as possible electrolytes of lithium-ion batteries in order to avoid several safety concerns present when organic electrolytes are used. Unfortunately, most tested IL pairs have been proven to be poor electrolytes in batteries [134,135]. One possible approach to improve specific IL ion transfer performance for the design of electrolyte-electrode couples in batteries includes a deeper understanding of thermodynamic and transport properties at the atomic-level. For that reason, computational simulations may be used not only to study these systems, but also to predict their properties, to narrow the wide variety of possible cation-anion combinations.
Some of us employed AMOEBA-IL to investigate the properties of an IL pair as a possible electrolyte candidate involving spyrocyclic pyrrolydinium ( was calculated for a range of temperatures between 300 to 500 K using the set parameters for one (1G) and three (3G) polarizable groups (Figure 2), as explained in Section 2. A decrease between 4%-8% can be observed in densities calculated considering intra-molecular polarization in the force field (3G) (Figure 12). Interestingly, the comparison between densities using 1G and 3G parameter sets shows a decrease of 12.7% and 16.1%, respectively, in a non linear fashion from 300 to 500 K. The 1G parameter set shows a density of 1.19 g/cm 3 and a volume of 297.45 Å 3 at 400 K, and 1.15 g/cm 3 (303.31 Å 3 ) at 450 K. On the other hand, for the 3G set of parameters, the densities and volumes at 400 K and 450 K were 1.14 g/cm 3 (324.70 Å 3 ) and 1.09 g/cm 3 (308.12 Å 3 ), respectively.
The mixture of these IL with 10% Li was also simulated. In general, the calculated densities for the mixture with both set of parameters are approximately 2%-4% higher than the density in the neat IL at the same range of temperatures ( Figure 12). These results are consistent with densities obtained experimentally for other neat and Li/IL mixtures [136][137][138][139]. Similar to the neat IL, a noticeable change in density for the 1G and 3G systems (5.4% and 3.1%, respectively) was observed in the mixture between 400 and 450 K. For the 1G set of parameters at 400 K, the density (volume) correspond to 1.21 g/cm 3 (275.06 Å 3 ), compared with 1.18 g/cm 3 (282.31 Å 3 ) at 450 K. The density(volume) for the 3G system at 400 K are 1.18 g/cm 3 (282.84 Å 3 ), whilst at 450 K the corresponding properties are 1.13 g/cm 3 (296.31 Å 3 ). The self-diffusion coefficients were also calculated for the same range of temperatures. The self-diffusion coefficients for the 1G system exhibits an increase of 5.6% in D ± from 300 K to 500 K, which is not considered a significant change compared with the D ± = 65.1% observed with the inclusion of intra-molecular polarization (3G) from 300 K to 500 K. These results suggest that a more accurate description of many body interactions speeds up the diffusion of the ions in the system. From 300-400 K, the 3G parameter set does not show a significant change in the self-diffusion coefficient (<1.5%), while an increment in D ± of 15.3% from 400 K to 500 K is observed. Conversely, the diffusion coefficient changes when intra-molecular polarization is neglected (1G) are not significant (1.7%). These results are consistent with the expected self-diffusion coefficients for smaller cations, showing that the bigger size of the [sPyr + ] cation results in slower diffusion. On the other hand, for the 10% Li + mixture, the ions diffuse faster using the 1G parameter set: for temperatures between 450 and 500 K there is an increment in the diffusion coefficient with respect the neat IL. The mixture with intra-molecular polarization (3G model), shows that self-diffusion coefficient decreased at the same range of temperatures compared to neat IL 3G model. The observed differences in calculated results between the 1G and 3G parameter sets for both the neat and Li mixtures are due the shortcomings of the 1G model in describing the change in charge density distribution caused by the changes in the internal structure of SPyr, specifically the bending of the cycles. Conversely, the 3G set enables the response of the polarization due to changes in the cation structure, providing a better representation of the inter-molecular interactions in the system. Radial Distribution Functions for anion-anion (B-B), anion-cation (B-N), and cation-cation (N-N) for both sets of parameters were calculated between 300 and 500 K as shown in Figure 13. The B-B RDF shows a small peak at 4 Å at 300, 350 and 400 K, which can be explained by the proximity of the [BF − 4 ] ions at those temperatures. Interestingly, this peak vanished at higher temperatures for the 3G model. The N-N RDF shows a plateau at 9.5 Å for the 300-400 K range, but it also vanished at higher temperatures. These differences support the hypothesis that the inter-ionic interactions are stronger with the inclusion of the intra-molecular polarization in this system. Additionally, the differences in RDFs at different temperatures, especially the loss of the peak for the B-B RDF at close distance, suggest a difference in the ordering of the ions at low temperatures compared with the high temperature systems. The above results suggest a significant change in thermodynamic and transport properties at temperatures between 400 K to 450 K. It should be noted that the computational results were obtained prior to the experimental synthesis and characterization of the neat [sPyr + ][BF 4 − ]. Subsequent to the observation in the significant changes in RDFs (and other properties) between 400 and 450 K, the neat IL pair was synthesized and differential scanning calorimetry (DSC) was performed. The inset in Figure 14 corresponds to the experimental DSC thermogram, which shows that this IL pair has a melting point of 448 K, enthalpy of fusion of 181 J/g, crystalization onset of 446 K and enthalpy of crystalization of 350 J/g. Thus, the calculated structural (and thermodynamic/transport) changes observed between 400 and 450 K are in very good agreement with the experimental results.

QM/MM Simulation of an Aniline Protection Reaction in H 2 O/[EMIm][BF 4 ]
AMOEBA-IL has also been applied to study chemical reactions in IL mixtures. In particular, a sample organic reaction of the N-tert-butoxycarbonylation of aniline in an IL mixture by polarizable QM/MM was recently performed [140]. To our knowledge, this is the first example of a QM/MM simulation of a reaction in ILs using multipolar/polarizable methods for the classical environment. This is an example of an important reaction to control functional groups in the synthesis of drug molecules, also known as N-t-Boc reactions [141]. Several reactions mechanisms, considering the effects of ionic liquids (water/1-ethyl, 3-methyl imidazolium/tetrafluoroborate ([EMIm][BF 4 ])/water mixture) as solvents via QM/MM simulations were investigated. The reaction mechanisms were characterized by performing minimum energy path (MEP) optimizations using two different chain-of-replica methods: the quadratic string method (QSM) and the nudged elastic band (NEB) [112,142].
Depending on which group is attacked, the reaction may take place via two different mechanisms (see Figure 15): The first corresponds to a step-wise mechanism, where the nucleophilic attack does not happen concommitantly with the formation of CO 2 . The second mechanism proceeds via an concerted path where CO 2 is formed at the same time as the nucleophilic attack occurs. These two possible mechanisms were studied using two different configurations for the system: configuration 1 was obtained by restraining the distance of one IL ion pair to the solute (aniline and BoC) in the QM/MM optimization of snapshots taken from MD simulations. This configuration includes 149 atoms in the QM region.
Configuration 2 was obtained based on QM/MM optimization of a series of snapshots from unrestrained MD simulations where the IL anions and cations are not close to the solvent. Configuration 2 includes 155 atoms in the QM subsystem. Both configurations were investigated in water/IL QM/MM systems. In the first mechanism, the nucleophilic attack in configuration 1 is followed by a proton transfer from the aniline to the Boc group, allowing the formation of tert-butyl carbamate and tert-butyl carboxylic acid as products ( Figure 16). For configuration 2, the CO 2 is formed concurrently with the carbamate. The energy barrier in this case is roughly 3 kcal/mol lower compared with configuration 1. Fluctuations in anisotropic field due to the presence of IL pairs stabilize a metastable intermediate (MI) produced by the nucleophilic attack, coupled with a second low-barrier TS corresponding to the proton transfer followed by the formation of the products (Figure 17). In the second mechanism, the minimum energy path for configuration 1 did not converge. On the other hand, configuration 2 shown in Figure 18 has a rate-limiting step barrier similar to the barriers of the configurations in mechanism 1, corresponding to the nucleophilic attack step. It also shows a (MI) previous to a second energy barrier that corresponds to the proton transfer of the aniline to the Boc groups that leads to the formation of tert-butanol, CO 2 and tert-butyl carbamate.
Analysis of the electronic wavefunctions via combined electron localization function (ELF) and non-covalent interaction index (ELF/NCI) were performed to obtain more insights related to the evolution of the intra-and intermolecular interactions in the system. The ELF/NCI analyses for the TS structures for both configurations for the rate-limiting step in mechanism 1 ( Figure 19) shows a disynaptic basin between the C1 and N43 which suggests the formation of a covalent bond between these two atoms in the TS. This means that in all cases, the TS structure corresponds to the mentioned nucleophilic attack and is part of a late TS. By contrast the TS structures for the second mechanism in both configurations did not show a disynaptic basin between N43 and C1. Nevertheless, a strong attractive interaction in the form of a ring can be observed in Figure 20b.
These results indicate a stabilization of an early TS provided by the solvent environment. Furthermore Figure 20c shows a bifurcated H-bond interaction between H45, O11 and O13 which explains the transfer of H49 to O11 in TS2. The NCI analyses did not show surfaces between H45 and N43 in TS2 which suggests a longer distance between the products.
Finally, MEPs for both mechanisms using configuration 2 were re-optimized neglecting the polarization of the classical region to investigate the effect of the polarization of the environment on the reaction path. The first mechanism resulted in an energy profile similar to the one obtained with explicit polarization. However, this MEP did not show the second TS corresponding to the proton transfer. The second mechanism showed an energy profile with some irregularities suggesting a rearrangement of environmental solvent molecules The energy barrier for this mechanism was almost twice the value of the MEP shown in Figure 18. The similarities observed between the energy profiles with and without explicit polarization presented in Figure 21 can be explained because of the use of the optimized path with polarization as initial guess for the MEP optimization without polarization, considering also that the optimized path corresponded to the potential energy surface rather than free energy surface.   These results show that the N-tert-butoxycarbonylation of aniline can take place via either concerted or sequential reaction mechanism, depending on the chosen configuration for the ionic pair. Overall, the rate-limiting step for this reaction corresponds to the nucleophilic attack of the aniline. The products for each mechanism depend on which Boc group is attacked for the aniline. ELF/NCI analyses suggest that for mechanism (1), where no CO 2 is formed, the rate-limiting step corresponds to a late Transition State (TS). On the other hand, for the second mechanism, where concomitant formation of CO 2 along the rate-limiting step occurs, an early TS structure is stabilized by the surrounding environment. The inclusion of explicit polarization of the MM environment demonstrated to provide a better representation of the the change density distribution of the QM subsystem, and that it has an important effect on the energetics and profile of the MEP.

Conclusions
AMOEBA-IL provides accurate results for several ionic liquid systems. Due to the highly charged nature of the ionic liquids, the inclusion of explicit polarization is essential in the description of different properties. Simulations of various IL systems using AMOEBA-IL have resulted in the reproduction of thermodynamic (heat of vaporization, density), transport (diffusion coefficients) and kinetic (water-exchange rates) properties, and prediction of others (phase-changes), that were later corroborated experimentally. However, challenges still remain, for example, in the case of the calculated exchange rate coefficients for the pure water and water/IL systems, the calculated values show deviation with respect to the experimental results, although the trends are well correlated with experiment. AMOEBA-IL's explicit intra-molecular polarization allowed a better representation of the distribution in charge density due to changes in the configuration of the [SPyr] cation, resulting in a better description of density and diffusion coefficients. Finally, the use of polarization in the MM-environment in QM/MM calculations has proved to be quite important to obtain reliable results in the calculation of energy barriers in systems using IL as solvents. AMOEBA-IL is a multipolar-polarizable potential for ILs that provides accurate description of inter-molecular interactions and bulk properties, and can be used for classical and hybrid QM/MM calculations. Future development of AMOEBA-IL will include expansion of the parameter library and applications to different systems using both classical and QM/MM simulations.