From molten salts to room temperature ionic liquids: Simulation studies on chloroaluminate systems
Mathieu Salanne, Leonardo J. A. Siqueira, Ari P. Seitsonen, Paul A. Madden, Barbara Kirchner
aa r X i v : . [ c ond - m a t . m t r l - s c i ] M a r From molten salts to room temperature ionic liquids: Simulation studies onchloroaluminate systems
Mathieu Salanne , Leonardo J. A. Siqueira , Ari P. Seitsonen , Paul A. Madden , and Barbara Kirchner UPMC Univ Paris 06, CNRS, ESPCI, UMR 7195, PECSA, F-75005 Paris,France. Fax: +33 1 4427 3228; Tel: +33 1 4427 3265; E-mail: [email protected] Laborat´orio de Materiais H´ıbridos, Universidade Federal de S˜ao Paulo, CEP 04024-002, S˜ao Paulo, SP, Brazil Physikalisch-Chemisches Institut, University of Zurich,Winterthurerstrasse 190, CH-8057 Zurich, Switzerland Department of Materials, University of Oxford, Parks Road, Oxford OX1 3PH, UK. and Wilhelm-Ostwald Institute of Physical and Theoretical Chemistry,University of Leipzig, Linn´estr. 2, D-04103 Leipzig, Germany
An interaction potential including chloride anion polarization effects, constructed from first-principles calculations, is used to examine the structure and transport properties of a series ofchloroaluminate melts. A particular emphasis was given to the study of the equimolar mixture ofaluminium chloride with 1-ethyl-3-methylimidazolium chloride, which forms a room temperatureionic liquid EMI + -AlCl − . The structure yielded by the classical simulations performed within theframework of the polarizable ion model is compared to the results obtained from entirely electronicstructure-based simulations: An excellent agreement between the two flavors of molecular dynamicsis observed. When changing the organic cation EMI + by an inorganic cation with a smaller ionicradius (Li + , Na + , K + ), the chloroaluminate speciation becomes more complex, with the formationof Al Cl − in small amounts. The calculated transport properties (diffusion coefficients, electricalconductivity and viscosity) of EMI + -AlCl − are in good agreement with experimental data. I. INTRODUCTION
When monovalent chloride salts are dissolved in liquidAlCl , ionic liquids with low melting points are formed.When the added cation is inorganic, the system becomesa molten salt, like NaAlCl with T m = 152 ◦ C . When itis organic, the system is likely to be liquid below 100 ◦ C,thus potentially becoming a room temperature ionic liq-uid (RTIL); for example T m (EMIAlCl ) = 7 ◦ C (whereEMI + stands for 1-ethyl-3-methylimidazolium cation).A number of investigations using molecular dynamicssimulations have been performed for ionic liquids directedat either structure or dynamical properties. Given themolecular nature of cations or anions forming ionic liq-uids, earlier molecular dynamics simulations were per-formed with pairwise additive potentials based on non-polarizable force fields such as OPLS or AMBER.
These simulations have been able to reproduce the struc-tural properties of ionic liquids, for instance showing thata somewhat unusual hydrogen bond can take place inimidazolium based ionic between the imidazolium ringproton and the anion species. This observation is inagreement with IR and NMR spectroscopy studies.
At intermediate spatial range (0.5 to 0.8 ˚A − ), a chain-length dependent pre-peak appears either in the staticstructure factor calculated from molecular dynamics tra-jectories or obtained by small wide angle X-ray scat-tering (SWAXS) and/or small angle neutron scattering(SANS). This feature has been assigned to structuralsegregation or inhomogeneity, i.e. the imidazolium ringsand the alkyl chain regions of cations respectively clusterto form polar and apolar domains in ionic liquids wherethe cations have long alkyl chains. Despite their ca-pability to reproduce the structure, the non-polarizable force fields often begin to fail when predicting dynami-cal properties, such as diffusion coefficient, electrical con-ductivity, and viscosity. In general the dynamical prop-erties calculated with non-polarizable models are slowerthan is observed in experiments: The diffusion coefficientand electrical conductivities are lower while the viscos-ity is greater than the measured values.
Note thatpairwise additive potentials are often parameterised to implicitly include some aspects of the polarization effects.For example, the simulated fluid may become less viscousif the magnitudes of the partial charges are reduced, but sometimes at the price of a poorer representationof other properties. These trends are particularly well-documented in the case of the inorganic systems. Disregarding the polarization method used, moleculardynamics simulations of molten salts in which polariza-tion effects are included have shown that the liquids be-come more mobile.
Moreover, molecular dynamicssimulations of ionic liquids including polarization haveshown better agreement of dynamical properties withexperiment, that is, the simulated liquids are either asmobile as the experimental sample or more mobile thanthe liquids simulated with non-polarizable model.
Among the studies using a polarizable model for ionicliquids, we would cite the recent and comprehensive workby Yan et al, where the authors carefully investigatedthe role-played by polarization in the ionic liquid EMI + NO − . The authors found that when polarization effectsare included, the short-range interactions are enhanceddue to the additional charge-dipole and dipole-dipole in-teractions, while long-range electrostatic interactions be-come more screened. Accordingly, the imidazolium ringproton - NO − hydrogen bond strength is increased uponswitching the polarization on. Another consequence ofthe stronger polar short-range interactions achieved withthe polarizable model for the structure is the enhance-ment of the tendency towards structural inhomogeneity,which acts by pulling the apolar chains away from thepolar center and reinforces the assembly of ethyl chaininto the apolar domains. The attenuation of long-rangeelectrostatic interaction is responsible for the increase ofdiffusion, conductivity and structural relaxation and thereduction of viscosity.In this work, we study a series of chloroaluminate liq-uids. A particular emphasis is given to the room tem-perature ionic liquid consisting of an equimolar mixtureof EMICl with AlCl . We employ two flavors of molecu-lar dynamics simulations, which differ by the method in-volved for the calculation of the interactions: In the first,they are obtained from a full density functional theory(DFT) calculation, and in the second an analytical forcefield (including polarization effects) is used. Both meth-ods have already been used successfully for the study ofpure AlCl . In a first step we detail how the param-eters of the interaction potential for EMICl-AlCl havebeen obtained. The structure and dynamics of this meltare then analyzed and compared to a series of moltensalts consisting of equimolar mixtures of AlCl with LiCl,NaCl and KCl. II. INTERACTION POTENTIALSA. Functional form
In recent work on the simulation of inorganic ionic sys-tems, interaction potentials that were used commonly in-cluded many-body effects such as polarization.
Themost generally retained analytical expression is V = X i As soon as a functional form has been chosen for theinteraction potential, all the atomic and atom pair pa-rameters have to be carefully determined. In all cases,the overall charges on the Cl − , Al , and EMI + are fixedas -1, +3, and +1 respectively. Our approach consistsin parameterizing the remaining terms on the basis offirst-principles electronic structure calculations. This isachieved by a “force-matching” method, which is gen-eralized by using the first-principles calculated induceddipoles in addition to the forces in the fitting procedure.Amongst molten chlorides, this method has already beenapplied successfully to the case of LiCl, NaCl and KCl. Here we extend this set of parameters in order to includethe Al and 1-ethyl-3-methylimidazolium cations. FIG. 1: Topology of the 1-ethyl-3-methylimidazolium molec-ular ion.TABLE I: Interaction potential parameters (partial charges,Lennard-Jones parameters) used for the EMI + cation. Notethat the ǫ i and σ i parameters were used for the cation-cationvan der Waals interactions only.Atom i ǫ i σ i q i (kJ mol − ) (˚A) ( e )N In the present work, we have chosen to describe theEMI + cation as a non-polarizable species for two reasons.Firstly, polarization effects are dominated by the chlorideanions (in systems like EMI-BF involving less polariz-able anions the situation is opposite, with the polariza-tion of the cation playing the main role ). Secondly,from a practical point of view, the computer time associ-ated with the minimization procedure of the polarizationterm increases substantially with the number of polar-izable species. For all the intramolecular/intermolecularinteractions involving the EMI + cation only, the param-eters were taken from the well-tested force field of P´aduaand co-workers. The topology of the molecule is shownon figure 1, and the partial charges and Lennard-Jones individual parameters are given in table I. The pair pa-rameters involved in equation 5 are obtained from theusual Lorentz-Berthelot mixing rules: ǫ ij = p ǫ i × ǫ j (6) σ ij = σ i + σ j FIG. 2: Results of force matching for one EMICl-AlCl con-figuration. The three components of the force acting on eachaluminium and chloride ions are compared (dots: DFT, line:PIM (dots: DFT, line: PIM). We have then applied the force-matching method tofit all the other interactions, which involve the chlo-ride anion on the one hand and either the aluminiumor the EMI + cation atoms on the other hand. In or-der to ensure that the force-matching really refines thecation-anion terms, we have held fixed the parametersfor the Cl-Cl interactions at the values previously de-termined for the pure alkali halides. To this end, weextracted several configurations from some trajectoriesobtained from Car-Parrinello molecular dynamics simu-lations of the liquids EMICl and EMICl-AlCl . Each ofthese configurations were then used again as an inputto a plane wave density functional theory (DFT) calcu-lation (the details of these calculations are given in thenext section). We first performed an energy minimiza-tion in order to find the ground-state electronic structure.The forces on each chloride and aluminium (if present)ions were extracted at this stage. Then we transformedthe Kohn-Sham orbitals into a set of maximally local-ized Wannier functions, from which we could determinethe induced dipole moment on the Cl − anions. The pa-rameters in the polarizable potentials are then optimizedby matching the dipole and forces from the potential tothe DFT calculated ones. Two distinct functions areminimized. The first one corresponds to the polarizationpart: χ D = 1 N Cl − N Cl − X i =1 | µ i DFT − µ i PIM | | µ i DFT | (8)where the subscript PIM stands for “Polarizable IonModel”. Since the chloride anion is the only polariz-able species, for a given atom type α the two short-rangedamping parameters b Cl α and c Cl α were optimized withthis function. The second function, χ F = 1 N Cl − N Cl − X i =1 | F i DFT − F i PIM | | F i DFT | (9)+ 1 N Al N Al3+ X i =1 | F i DFT − F i PIM | | F i DFT | , is used to parameterize the repulsion term of the inter-action potential, i.e. for determining the various B ij and a ij values in equation 2. A well known limitation of thePBE functional that we used in our DFT calculation isthat it does not take into account properly the disper-sion interactions, so that the C ij and C ij cannot befitted by this procedure. Instead, we have used our pre-viously determined Cl − -Cl − dispersion coefficients forcalculating the various C ijn ( n = 6,8) by using the Slater-Kirkwood relationship.The quality of the representation of the forces is il-lustrated in figure 2 where we compare the three com-ponents of the forces acting on the aluminium (32 firstatoms) and chloride ions calculated with our potentialwith the DFT calculated data, for one of the EMICl-AlCl configurations. The agreement is seen to be uni-formly good. For both RTIL systems, we obtained valuesof 0.13 for χ F , which is similar to our previous work onmolten LiCl, NaCl and KCl. Concerning the dipoles, thevalues of χ D are respectively 0.39 and 0.03 on average forthe EMICl and EMICl-AlCl configurations. The pooreragreement in the first case is due to the fact that theelectric field felt by the Cl − anions, which is only dueto the EMI + cation, is only partially reproduced by thepartial charges of the atoms. When aluminium cationsare added to the melt, this electric field becomes domi-nated by the trivalent charges they carry, and this effectbeing much more easily captured by our model. The ob-servation suggests that further refinement of the chargesof the EMI + model would be a good target to furtherimprove the predictive power of the simulations.The parameters which were determined from the fit-ting procedure are summarized in tables II and III. Nodistinction has been made between the various H i , C i and N i atom types because it did not lead to any im-provement of the fit quality. TABLE II: Parameters for the pairwise additive terms of theinteraction potential. b ij = b ij = 3.21 ˚A − for all pairs.Atom pair B ij × − a ij C ij × − C ij × − i − j (J mol − ) (˚A − ) (J mol − ˚A ) (J mol − ˚A )Cl-Cl 722.0 3.396 8.07 4.52Cl-Al 95.8 2.955 0.0 0.0Cl-Li 154.5 3.685 0.0 0.0Cl-Na 177.2 3.262 2.16 2.70Cl-K 547.5 3.458 1.52 1.91Cl-C i i i − , Na + and K + ions were respectively set to 2.96, 0.13 and 0.70 ˚A .Atom pair b ij c ij c ji i − j (˚A − ) (–) (–)Cl-Al 3.802 2.953 –Cl-Li 3.517 2.079 –Cl-Na 3.326 3.000 0.697Cl-K 3.084 3.000 0.917Cl-H i III. SIMULATION DETAILSA. Electronic structure-based molecular dynamics For the ab initio-molecular dynamics simulations weused density functional theory with the generalised gra-dient approximation of Perdew, Burke and Ernzerhof,PBE as the exchange-correlation term in the Kohn-Sham equations, and we replaced the action of the coreelectrons on the valence orbitals with norm-conservingpseudo potentials of the Troullier-Martins type. Weexpanded the wave functions in plane waves up to thecut-off energy of 70 Ry. We sampled the Brillouin zoneat the Γ point, employing periodic boundary conditions.We performed the simulations in the NVT ensem-ble, employing a Nos´e-Hoover thermostat at a targettemperature of 300 K and a characteristic frequency of595 cm − , a stretching mode of the AlCl molecules.We propagated the velocity Verlet equations of motionwith a time step of 5 a.t.u. = 0.121 fs, and the ficti-tious mass in the Car-Parrinello molecular dynamics for the electrons was chosen as 700 a.u.; we note thatthe fictitious dynamics increases the temperature seenby the ions somewhat. A cubic simulation cell with aedge length of 22.577 ˚A containing both 32 anionic andcationic molecules, equalling to the experimental densityof 1.293 g/cm . The length of the trajectory was 20 ps.The CPMD simulation code was used. B. Polarizable ion model-based molecular dynamics For the inorganic molten salts (LiCl-AlCl , NaCl-AlCl and KCl-AlCl ) the simulations were performedin the NVT ensemble, employing a Nos´e-Hoover ther-mostat at a target temperature of 573 K, with a timestep of 0.5 fs, for a total simulation time of 3 ns. ForEMICl-AlCl , several simulations were performed at thetemperatures of 298, 433, 473, 523 and 573 K, with a timestep of 1 fs. Note that the higher temperatures employedmay be above the thermal decomposition temperature ofthe ionic liquid (no data could be found in the littera-ture); these simulations were performed in order to makea fair comparison of the structure and dynamics with theinorganic molten salts. The total simulation time was of1.5 ns for each temperature except 298 K, for which itwas of 3 ns. The FIST module of CP2K simulation pack-age was used. In all cases, cubic simulation cells wereemployed. The total number of atoms and the cell edgelength, which were chosen according to the experimentaldensities (an extrapolation of the room temperature datawas made for EMICl-AlCl when necessary), , are pro-vided in table IV. The cut-off distance for the real spaceof the Ewald sum and the short range potential was setto 13 ˚A for EMICl-AlCl and to half the size of the sim-ulation cell otherwise. TABLE IV: Number of atoms and cell edge length ( L ) em-ployed in the simulations ( M + = Li + , Na + , K + or EMI + ).System T (K) N M + N Al N Cl − L (˚A)LiCl-AlCl 573 120 120 480 28.426NaCl-AlCl 573 120 120 480 28.790KCl-AlCl 573 120 120 480 29.728EMICl-AlCl 298 200 200 800 41.576433 200 200 800 42.807473 200 200 800 43.200523 200 200 800 43.713573 200 200 800 44.251 IV. RESULTS AND DISCUSSIONA. Chloroaluminate speciation Unlike previous classical MD simulations ofchloroaluminate-based ionic liquids, in which theAlCl − ion was treated as a molecular species, anadvantage of the two simulation methods involved in thepresent study is that they both treat the chloride andaluminium ions as independent species. This is doneintrinsically in electronic structure based simulationson one hand, and by construction of the polarizableion model (PIM) on the other. This capability hasbeen used, for example, in pure AlCl liquid, in whichstructure is characterized by the formation of AlCl − tetrahedra sharing an edge to form Al Cl dimers orlarger clusters. MD simulation studies either employingsimilar interaction potentials as in the present study or based on electronic structure calculations wereable to confirm quantitatively this picture, in agreementwith X-ray diffraction experiments. The formation oflarge chloroaluminate species was also investigated insimulations of a system consisiting of one EMICl pairdissolved in a AlCl solvent. In equimolar mixtures of AlCl with monovalent chlo-ride M Cl salts, the main “chemical reaction” involvingchloroaluminate species which is expected to occur is: AlCl − + AlCl − → Al Cl − + Cl − (10)Here in the case of the equimolar EMICl-AlCl thetwo simulation methods agree on the formation of iso-lated AlCl − anions only, in agreement with experimen-tal results. Small differences are observed in the firstneighbour distances, which are shorter in the simulationemploying the polarizable model: The first peaks of theAl-Cl radial distribution functions respectively appear at2.18 ˚A and 2.22 ˚A for PIM-MD and CPMD, while the cor-responding Cl-Cl functions display some first peak max-ima at 3.53 ˚A and 3.61 ˚A. These differences may be at-tributed to the choice of keeping the Cl − -Cl − interactionpotential parameters similar to the values which had beenpreviously been determined for pure LiCl, NaCl and KClmolten salts. Another source of difference could be thepresence of attractive dispersion effects in the classicalMD only. FIG. 3: Al -Cl − radial distribution functions for the M Cl-AlCl mixtures at 573 K. Insert: Enlargment of the first min-imum zone. When M + is an inorganic cation (Li + , Na + , K + ),we observe a weakening of the Al-Cl links. This maybe quantified from the corresponding radial distributionfunctions (RDF), for which we observe a decrease of theintensity of the first peak, which switches from a valueof 47 in EMICl-AlCl to one of 18 in LiCl-AlCl (notshown). Some Cl − ions begin to jump from one alu-minium coordination sphere to another. The signatureof these exchanges can also be tracked in the RDFs whichare given on figure 3: The minimum occurring after thefirst maximum no longer corresponds to a value of zero asin the case in EMICl-AlCl . The value of this minimumalso increases when the ionic radius of the M + cation de-creases in the series. The position of the second peak ofthe RDF is slightly shifted toward larger distances whenpassing from Li + to Na + and then K + , but this distanceincreases by as much as ≈ + toEMI + (note that the peak also becomes much wider, witha pronounced shoulder). The absence of exchange of Cl − between AlCl − coordination tetrahedra in EMICl-AlCl can therefore easily be explained by the longer distancethat has to be crossed during a jump. FIG. 4: Al -Al radial distribution functions for the M Cl-AlCl mixtures at 573 K. Going back to the chemical reaction defined in equa-tion 10, we observe some signature of the formation ofAl Cl − ions in our systems by examining the Al-Al RDFs(shown on figure 4). A first peak with a small intensityis obtained in all the systems except EMICl-AlCl . Thecorresponding distance (ranging from 3.5˚A in LiCl-AlCl to 3.7 ˚A in KCl-AlCl ) is smaller than the length of twoAl-Cl bonds, indicating pairs of Al ions linked by com-mon Cl − anions. To quantify the proportion of eachchloroaluminate species present in the melt, we have fol-lowed the same procedure as in our previous work onLiF-BeF mixtures, where we could provide a quanti-tative picture of the speciation in the melt for a widerange of composition. Here we have chosen to definean Al-Cl-Al bond when two Al-Cl and the correspond-ing Al-Al distances are shorter than the correspondingRDFs minima. By extending the analysis of the linkagesbetween ions, we can calculate the proportion of any dif-ferent polynuclear ionic species present in the melt. Themain species that we observed apart from the isolatedAlCl − anions were the Al Cl − and Al Cl − species; their proportions are provided in table V for a temperature of573 K. Note that these numbers have to be taken cau-tiously: Due to the finite size of the simulation cells andto the relatively low number of chloride anions jumpingor bridging events, some refinement may be necessary.We also observed some small proportions of transientchemical species such as AlCl − or Al Cl − in negligi-ble amounts. TABLE V: Proportions (given in percentage of aluminiumions) of the various chloroaluminate species at 573 K.System [AlCl − ] [Al Cl − ] [Al Cl − ]LiCl-AlCl B. Effect of the M + cation nature on theintermolecular structure FIG. 5: AlCl − - M + radial distribution functions at 573 K. Inthe case where M + is the EMI + cation, the position of thecenter of mass of the molecule was used to calculate the func-tion. Small arrows indicate the position of the correspondingAl -Al RDF first maximum (occulting the first peak dueto the formation of Al Cl − ions); the Na + one is not visiblebecause it is confounded with the K + one. From the analysis of speciation in the chloroaluminatespecies we can try to describe the equimolar mixtures of M Cl with AlCl as generic M X molten salts (like e.g.alkali halides), where M + is the monovalent cation and X − is the AlCl − anion. In simple M X molten salts, thestructure around a given ion consists in alternating layersof opposite charges. As a result, the RDFs are character-ized by strong oscillations, with the like-like RDFs max-ima (minima) which are located at the same positionsas the cation-anion RDF minimum (maximum). In thepresent work, we can compare the M + -Al RDF to theAl -Al one; the Al was chosen because it gives theAlCl − ion center of mass; for the EMI + cation the M + -Al was calculated with respect to the imidazolium ioncenter of mass position. These functions are shown onfigures 4 and 5. On the latter we also report the positionof the Al -Al RDF first maximum (omitting the smallpeak due to the formation of Al Cl − ) with an arrow onthe abcissa in order to facilitate the comparison with theposition of the first minimum in the cation-anion term.We observe that these arrows are not located at the sameposition as the corresponding minima, but the differenceremains small. This means that the structure deviatesslightly from that of a simple M X molten salt. Thisconclusion still holds for EMICl-AlCl at lower temper-atures: The main features of the M + -Al RDF remainunchanged when passing from 573 K to 298 K.From this figure we can observe somewhat differentbehaviour depending on the size of the M + cation. In thecase of the LiCl-AlCl melt, the M + -Al RDF minimumis located at a shorter distance than the correspondingAl -Al RDF first maximum; this is due to the smallsize of the Li + cation, which allows it to “penetrate”slightly into the AlCl − anion van der Waals sphere. Thedifference becomes smaller for the NaCl-AlCl and KCl-AlCl melts for which the extrema of the two functionsalmost coincide, and the situation is reversed for EMICl-AlCl where closer packing of the anions is allowed. The M + -Al RDF first peak is also broader in the case ofEMICl-AlCl compared to the alkali cations containingsystems, but the difference is not dramatic. The twolatter observations suggest that in the present system theEMI + cation behaves mostly like a monoatomic cationwith a very large ionic radius. In the next section wewill investigate its coordination structure on a more localscale, in order to understand if the charge delocalizationas well as the presence of different atom types induce aparticular ordering in the first solvation shell. C. EMICl-AlCl equimolar mixture In order to characterize further the intermolecularstructure of the EMICl-AlCl system, the radial distribu-tion functions of the chloride anions and the atoms fromthe EMI + cations have been calculated. They are plot-ted on figures 6 (hydrogen atoms) and 7 (carbon atoms).The CPMD and PIM results obtained at room tempera-ture are compared; an excellent agreement is observed forall the atoms. The small differences observed are morelikely due to the poorer sampling in the CPMD simula-tion rather than to a deficiency of the classical interactionpotential. In a previous attempt to obtain classical forcefield parameters through a force-matching procedure sim-ilar to ours on the ionic liquid dimethylimidazolium chlo-ride (DMICl), Youngs et al. have been less successful inreproducing the structure obtained from CPMD. Al-though they were able to reproduce the peak positions FIG. 6: Radial distribution functions of the chloride anionsand the protons from the EMI + cations. Top: CPMD results;bottom: PIM results. Left: Protons from the imidazoliumring. Right: Protons from the alkyl chains.FIG. 7: Radial distribution functions of the chloride anionsand the carbon atoms from the EMI + cations. Top: CPMDresults; bottom: PIM results. Left: Carbon atoms from theimidazolium ring. Right: Carbon atoms from the alkyl chains. better than with previous classical MD potentials, theintensities of the first peaks where largely overestimated(especially in the case of the ring protons). The over-structuring obtained by these authors could partly becancelled by multiplying all the ǫ i parameters by a fac-tor of 2, at the price of a worse reproduction of the ini-tial set of DFT forces (therefore leading to higher χ F values). The better agreement observed here can be at-tributed to two factors. The main one is the use of amore complicated functional form for the analytical in-teraction potential, or in other words to the inclusion ofanion polarization effects. The second is the fact that wefocused on the chloride (and aluminium) ion forces in ourforce-fitting process; in Youngs et al. study all the pa-rameters, including the intramolecular ones, were fittedduring the procedure. The intramolecular contributionmay have somewhat “masked” the intermolecular ones,thus leading to a smaller precision for the latter. It isworth noticing that the better agreement observed herewas not obtained by fitting more parameters. Due to thechoice of using the set of partial charges of Padua and co-workers for the EMI + ion atoms, only twelve parametersinvolving the Cl − and EMI + species had to be fitted inthe parameterization process: The eight Cl-C i ,N i ,H i ,Alrepulsion parameters ( B ij and a ij ) and the four param-eters b ij and c ij involved in the damping of the chlo-ride ions induced dipoles by the hydrogen and aluminiumatoms. FIG. 8: The electrostatic potential mapped onto the electrondensity with a surface values of 0.067 e − .˚A − . The colourscale ranges from -0.1 (red) to +0.1 (blue) atomic units. Left:AlCl − , right: EMI + . A striking feature of the radial distribution functionsinvolving the chloride ions and the EMI + atoms is thatthe most intense peaks are observed for the carbon atoms.The peak is more pronounced for the C and C , whileits intensity is almost the same for the other carbons.This result contrasts sharply with some other RTILs,like the above-mentioned DMICl for example, for whichthe anions tend to stay in the vicinity of the imida-zolium ring protons. On the contrary, it resembles tothe case of the 1- n -butyl-3-methylimidazolium hexafluo-rophosphate, for which the hydrogen bonds were foundto be weaker than expected, as indicated by their shortlifetimes. In the present system the situation is dueto the huge electrostatic interaction between Cl − anionsand the (Lewis) acidic Al cation. To illustrate this theelectrostatic potential mapped onto the isosurface of theelectron density of the two individual ions is shown infigure 8. It provides insight into the charge distributionof each ion. For the AlCl − we recognize in the red color(low electrostatic potential) the consequence of the nega- tive charge that is associated with this ion. The negativecharge is distributed all over the chlorine atoms, witha higher negative charge on the inner side (the sphere-like surface around the chloride is coloured in dark red onthe interior and in lighter red at the exterior of the AlCl − ion). Towards the center, close to the aluminium atom,we find a decreased negative charge shown as the bluecolor in the left panel of Fig. 6. The opposite is the casefor the EMI + . Here the positive charge leads to the bluecolor (high electrostatic potential) around this ion. Uponcloser inspection we find that the ring protons show thesame blue color as most of the molecule. A slight decreaseof the charges can be found in the methyl group protonsand a stronger decrease (green color) can be found at theterminal ethyl protons. This is in accordance with theobservation from the radial pair distribution functionsand with chemical intuition that these protons are lessacidic than the other protons of the cation.In the PIM simulations, the chloride anions thereforehave their induced dipoles oriented along the Al-Cl axis,which hinders the formation of hydrogen bonds. If welook in further detail at the chloride-proton radial dis-tribution functions, we see that the H -Cl − one showsthe most pronounced peak, associated to the shortestseparation. It is therefore clear that this proton is themost acidic coordination site. The functions involvingthe other two protons of the imidazolium ring (H andH ) also display some peaks at slightly larger distances.However, considering the protons from the ethyl and themethyl group, small pronounced peaks can also be ob-served. The methyl group hydrogen (H ) function firstmaximum has the same height as the H one; a similarobservation is made for the terminal hydrogen atoms ofthe ethyl group (H ) compared to the ring H protons.The arrangement of the chloride anions around the im-idazolium cations therefore seems to correspond more toan involved network including the alkyl chains as wellas the ring of the imidazolium cation, rather than to in-dividual ion pairs, confirming the results obtained formany imidazolium-based ionic liquids. In order toconfirm this picture, the three-dimensional density mapwas calculated for the chloride anions located in the firstsolvation shell of the EMI + . Views of these maps perpen-dicular and parallel to the imidazolium ring are providedon figure 9. Red, green and blue regions correspond re-spectively to high, intermediate and low probability offinding an anion around a given cation. By high proba-bility we mean higher than 80 % of the maximum densityof occurences in these figures, and low probability meansin between 40 % and 60 % of the maximum density ofoccurences (probability lower than 40 % is not shown).The important features of the radial distribution func-tions, i.e. a preference for some regions close to the C and C carbons and an absence of important correlationswith the protons, are recovered. In general, localizationeffects are much less important than in the DMICl case,and the formation of hydrogen bond is not the main influ-ence governing the structure of EMICl-AlCl . It is worth FIG. 9: Probability density maps of location of nearest-neighbor anions around the imidazolium ring in equimolarEMICl-AlCl mixture at 298 K. The left panel shows theview above the imidazolium ring plane, and the bottom panelshows the view on the ring plane. Red, green and blue regionscorrespond respectively to high, intermediate and low prob-ability of finding an anion around a given cation. By highprobability we mean higher than 80 % of the maximum den-sity of occurences, and low probability means in between 40 %and 60 % of the maximum density of occurences (probabilitylower than 40 % is not shown). mentioning that in a study involving imidazolium cationswith a series of halide anions of different size, the analy-sis of density maps showed that the larger the anion thelower is the probability of the existence of hydrogen bondbetween the anion and the cation protons. The presentresults agree with this view when considering that therelevant anion in our system is the full AlCl − entity.The structural analysis of the equimolar EMICl-AlCl mixture shows that it can almost be viewed as a classical M X molten salt, where M + = EMI + and X − = AlCl − ,with a structure characterized by a competition betweenclose packing and charge ordering effects occuring be-tween the two species. Small deviations from the M X system are observed due to the existence of weak spa-tial correlations arising from the formation of hydrogenbonds between the ring protons and the chloride anions. D. Transport In the recent simulation work on the room tempera-ture ionic liquids, it was shown that polarization effectsare important for reproducing accurately the transport properties. Here from the Einstein relation D α = lim t →∞ t h| δ r i ( t ) | i (11)where δ r i ( t ) is the displacement of a typical ions ofspecies α in time t , we obtained a diffusion coefficientof 1.04 × − m s − for the cation at 298 K. Thiscompares very well with the experimental value mea-sured by pulse-field gradient NMR spectroscopy, whichis 0.95 × − m s − at 293 K. For the anion we ob-tained a diffusion coefficient of 4.45 × − m s − at298 K. FIG. 10: Viscosities (top) and electrical conductivities (bot-tom) of EMICl-AlCl in the form of an Arrhenius plot. Thesimulated values, although obtained in a different tempera-ture regime, compare well with the experimental ones. From the mean squared displacement of the overallcharge density it is also possible to calculate the elec-trical conductivity of a melt according to the followingformula: σ = e k B T V lim t →∞ t h| X i q i δ r i ( t ) | i , (12)where e is the elementary charge and q i the magnitudeof the (partial or formal) charge on atom i . Notice thatthis expression involves the correlations between the dis-placements of different ions, which is essential in this caseas the motion of the Al and Cl − ions within a given0AlCl − is strongly correlated. The viscosity is calculatedby integration of the correlation function of the stresstensor: η = 1 k B T V Z ∞ h σ αβ (0) σ αβ ( τ ) i d τ, (13)where σ αβ is one of the components of the stress tensor.The two latter quantities, which are collective transportcoefficients (compared to the diffusion coefficients, whichare individual ), require much longer simulation times toget converged values. Such quantities may therefore onlybe calculated with the analytic force-field model, they arewell out of reach of full ab initio simulations. In the caseof EMICl-AlCl , we could extract their values for all thetemperatures except at 298 K. They are shown on figure10 in the form of an Arrhenius plot. On the same figurewe also show the experimental values; although thosewere obtained at lower temperatures, we observe a goodmatch between the two sets of data, which confirms thatour polarizable ion model yields the correct dynamics ofthe system. TABLE VI: Diffusion coefficients of the equimolar M Cl-AlCl systems at 573 K ( M + = Li + , Na + , K + or EMI + ). All valuesare in 10 − m .s − System D M + D Al D Cl − LiCl-AlCl At the temperature of 573 K, which is above the melt-ing point of the inorganic chloroaluminates, it is also pos-sible to compare the diffusion coefficients of the variousions in all the systems. They are reported in table VI.For all the species, we observe an increase of the diffu-sion coefficient with the M + cation size. Note that thenumber density also changes when switching from onesystem to another, which is likely to be the main causefor this evolution. Still it is possible to come to the con-clusions that in all systems the chloroaluminate aniondiffuses more slowly than the corresponding M + cation.We can also conclude that the organic cation-based ionicliquid, EMICl-AlCl , fits into the same pattern of be-haviour as the inorganic cations-based compounds in re-gard to the variation of transport properties with cationsize. V. CONCLUSION We have shown in this study that the polarizable ionmodel developed for simple molten salts can success-fully be transfered to chloroaluminate ionic liquids. The interaction potential was parameterized purely from first-principles through a force and dipole-matching proce-dure, and a common set of parameters could be ob-tained for the mixtures of AlCl with several inorganicsalts (LiCl, NaCl and KCl) as well as with the organic1-ethyl-3-methylimidazolium chloride. An advantage ofour model is that it allows us to quantify the chloroalu-minate speciation due to the choice of treating Al andCl − as independent ionic species.A particular emphasis was given to the study of theequimolar mixture EMICl-AlCl , for which the structureyielded by the classical simulations performed within theframework of the polarizable ion model is compared tothe results obtained from entirely electronic structure-based simulations: An excellent agreement between thetwo flavors of molecular dynamics is obtained, which en-ables us to conclude that the structure of EMICl-AlCl can almost be viewed as the one of a classical M X moltensalt, where M + = EMI + and X − = AlCl − . The structureis characterized by a competition between close pack-ing and charge ordering effects occuring between the twospecies. Small deviations from the M X system are ob-served due to the existence of weak spatial correlationsarising from the formation of hydrogen bonds betweenthe ring protons and the chloride anions. When decreas-ing the cation ionic radius (EMI + → K + → Na + → Li + )small proportions of Al Cl − due to the reactionAlCl − + AlCl − → Al Cl − + Cl − (14)start to be observed (Al Cl − is also observed in equimo-lar LiCl-AlCl ). An interesting extension of this workwould be to study the same series of systems undernon-stoechiometric conditions. The influence of the M + cation chemical nature in the Lewis acidity of themelt can also be quantified using a recently developedmethod. In future work we will try to perform suchcalculations in order to relate it to the chloroaluminatespeciation determined in the present study.The only many-body effect included in the interac-tion potential for EMICl-AlCl is the anion polarization,which enabled us to perform long enough runs to getthe dynamic properties of the system, which are outof reach of first-principles simulations. The calculatedtransport properties were compared to the experimentaldata, showing again a good agreement for both the in-dividual (diffusion coefficients) and collective (electricalconductivities and viscosities) dynamics of the system. Acknowledgments MS is grateful to the FAPESP and the Federal Univer-sity of S˜ao Paulo for a travel grant which enabled this col-laboration. LJAS also thanks FAPESP, proc. 08/08670-7.1 H. A. Oye. Metall. Mater. Trans. B , 31B:641–650, 2000. Jr. A. A. Fannin, D. A. Floreani, L. A. King, J. S. Landers,B. J. Piersma, D. J. Stech, R. L. Vaughn, J. S. Wilkes, andJ. L. Williams. J. Phys. Chem. , 88(12):2614–2621, 1984. E. J. Maginn. Acc. Chem. Res. , 40:1200–1207, 2007. S. L. Price, C. G. Hanke, and R. M. Lynden-Bell. Mol.Phys. , 99:801–809, 2001. J. de Andrade, E. S. B¨oes, and H. Stassen. J. Phys. Chem.B , 106(14):3546–3548, 2002. C. J. Margulis, H. A. Stern, and B. J. Berne. J. Phys.Chem. B , 106:12017, 2002. T. I. Morrow and E. J. Maginn. J. Phys. Chem. B ,106:12161–12164, 2002. S. M. Urahata and M. C. C. Ribeiro. J. Chem. Phys. ,120(4):1855–1863, 2004. J. N. Canongia Lopes, J. Deschamps, and A. A. H. P´adua. J. Phys. Chem. B , 108:2038, 2004. M. Salanne, C. Simon, and P. Turq. J. Phys. Chem. B ,110(8):3504–3510, 2006. L. J. A. Siqueira and M. C. C. Ribeiro. J. Phys. Chem. B ,111:11776, 2007. B. L. Bhargava and S. Balasubramanian. J. Chem. Phys. ,127:114510, 2007. A. Wulf, K. Fumino, D. Michalik, and R. Ludwig. ChemPhysChem , 8:2265, 2007. K. Fumino, A. Wulf, and R. Ludwig. Angew. Chem., Int.Ed. , 47:3830, 2008. K. Fumino, A. Wulf, and R. Ludwig. Angew. Chem., Int.Ed. , 47:8731, 2008. S. M. Urahata and M. C. C. Ribeiro. J. Chem. Phys. ,124:074513, 2006. H. V. R. Annapuredy, H. K. Kashyap, P. De Biase, andC. J. Margulis. J. Phys. Chem. B , 114(50):16838–16846,2010. A. Triolo, O. Russina, B. Fazio, R. Triolo, and E. Di Cola. Chem. Phys. Lett. , 457:362, 2008. C. Hardacre, J. D. Holbrey, C. L. Mulland, T. G. A.Youngs, and D. T. Bowron. J. Chem. Phys. , 133:074510,2010. O. Russina, A. Triolo, L. Gontrani, R. Caminiti, D. Xiao,L. G. Hines, R. A. Bartsch, E. L. Quitevis, N. Plechkova,and K. R. Seddon. J. Phys.: Condens. Matter , 21:424121,2009. B. L. Bhargava, R. Devane, M. L. Klein, and S. Balasub-ramanian. Soft Matter , 3:1395–1400, 2007. R. M. Lynden-Bell and T. G. A. Youngs. J. Phys.: Con-dens. Matter , 21(42):424120, 2009. C. Cadena, Q. Zhao, R. Q. Snurr, and E. J. Maginn. J.Phys. Chem. B , 110:2821, 2006. S. Tsuzuki, W. Shinoba, H. Saito, M. Mikami, H. Tokuda,and M. Watanabe. J. Phys. Chem. B , 113:10641, 2009. T. Koddermann, D. Paschek, and R. Ludwig. ChemPhysChem , 8:2464, 2007. N.-T. Van-Oanh, C. Houriez, and B. Rousseau. Phys.Chem. Chem. Phys. , 12:930–936, 2010. F. M¨uller-Plathe and W. F. van Gunsteren. J. Phys.Chem. , 103:4745–4756, 1995. B. L. Bhargava and S. Balasubramanian. J. Chem. Phys. ,123:144505, 2005. D. Marrocchelli, M. Salanne, and P. A. Madden. J. Phys.:Condens. Matter , 22(15):152102, 2010. P. A. Madden and M. Wilson. J. Phys.: Condens. Matter ,12:A95, 2000. L. J. A. Siqueira, S. Urahata, and M. C. C. Ribeiro. J.Chem. Phys. , 119:8002, 2003. T. Yan, C. J. Burnham, M. G. Del P´opolo, and G. A. Voth. J. Phys. Chem. B , 108(32):11877–11881, 2004. C. Schroder and O. Steinhauser. J. Chem. Phys. ,133:154511, 2010. D. Bedrov, O. Borodin, Z. Li, and G.D. Smith. J. Phys.Chem. B , 114:4984–4997, 2010. T. Yan, Y. Wang, and C. Knox. J. Phys. Chem. B ,114:6905, 2010. T. Yan, Y. Wang, and C. Knox. J. Phys. Chem. B ,114:6886, 2010. F. Hutchinson, M. Wilson, and P. A. Madden. Mol. Phys. ,99(10):811–824, 2001. B. Kirchner, A.P. Seitsonen, and J. Hutter. J. Phys. Chem.B , 110(23):11475–11480, 2006. P. Tangney and S. Scandolo. J. Chem. Phys. ,117(19):8898–8904, 2002. C. Merlet, P. A. Madden, and M. Salanne. Phys. Chem.Chem. Phys. , 12:14109–14114, 2010. M. Salanne, D. Marrocchelli, C. Merlet, N. Ohtori, andP. A. Madden. J. Phys.: Condens. Matter , 23(10):102101,2011. V. Bitrian, O. Alcaraz, and J. Trullas. J. Chem. Phys. ,134(4):044501, 2011. L. Giacomazzi, P. Carrez, S. Scandolo, and P. Cordier. Phys. Rev. B , 83:014110, 2011. R. J. Heaton, R. Brookes, P. A. Madden, M. Salanne,C. Simon, and P. Turq. J. Phys. Chem. B , 110(23):11454–11460, 2006. K. T. Tang and J. P. Toennies. J. Chem. Phys. ,80(8):3726–3741, 1984. N. Ohtori, M. Salanne, and P. A. Madden. J. Chem. Phys. ,130(10):104507, 2009. O. Borodin. J. Phys. Chem. B , 113(33):11463–11478, 2009. N. Marzari and D. Vanderbilt. Phys. Rev. B , 56(20):12847–12865, 1997. A. Aguado, L. Bernasconi, S. Jahn, and P. A. Madden. Faraday Discuss. , 124:171–184, 2003. S. Zahn and B. Kirchner. J. Phys. Chem. A , 112:8430–8435, 2008. J. P. Perdew, K. Burke, and M. Ernzerhof. Phys. Rev.Lett. , 77:3865–3868, 1996. N. Troullier and J. L. Martins. Phys. Rev. B , 43:001993,1991. N. Troullier and J. L. Martins. Phys. Rev. B , 43:008861,1991. R. Car and M. Parrinello. Phys. Rev. Lett. , 55(22):2471–2474, 1985. P. Tangney. J. Chem. Phys. , 124(4):044111, 2006. The CPMD consortium. CPMD version 3.13.2. CP2K developers group. G. J. Janz. J. Phys. Chem. Ref. Data , 17(2):1–309, 1988. J. de Andrade, E. S. B¨oes, and H. Stassen. J. Phys. Chem.B , 106(51):13344–13351, 2002. F. Hutchinson, M. K. Walters, A. J. Rowley, and P. A.Madden. J. Chem. Phys. , 110(12):5821–5830, 1999. A. L. L. East and J. Hafner. J. Phys. Chem. B , 111:5316–5321, 2007. R. L. Harris, R. E. Wood, and H. L. Ritter. J. Am. Chem.Soc. , 73:3151–3155, 1951. B. Kirchner and A.P. Seitsonen. Inorg. Chem. , 46(7):2751–2754, 2007. Jr. A. A. Fannin, L. A. King, J. A. Levisky, and J. S.Wilkes. J. Phys. Chem. , 88:2609–2614, 1984. M. Lipsztajn and R. A. Osteryoung. J. Electrochem. Soc. ,132(5):1126–1130, 1985. C. L. Hussey, T. B. Scheffler, J. S. Wilkes, and Jr.A. A. Fannin. J. Electrochem. Soc. , 133(7):1389–1391,1986. M. Salanne, C. Simon, P. Turq, R. J. Heaton, and P. A.Madden. J. Phys. Chem. B , 110(23):11461–11467, 2006. M. Salanne, C. Simon, P. Turq, and P. A. Madden. J.Phys. Chem. B , 111(18):4678–4684, 2007. T. G. A. Youngs, M. G. Del P´opolo, and J. Kohanoff. J.Phys. Chem. B , 110:5697–5707, 2006. M. G. Del P´opolo, R. M. Lynden-Bell, and J. Kohanoff. J.Phys. Chem. B , 109:5895–5902, 2005. W. Zhao, F. Leroy, B. Heggen, S. Zahn, B. Kirchner,S. Balasubramanian, and F. Muller-Plathe. J. Am. Chem.Soc. , 131(43):15825–15833, 2009. M. Kohagen, M. Brehm, J. Thar, W. Zhao, F. M¨uller-Plathe, and B. Kirchner. J. Phys. Chem. B , 115(4):693–702, 2011. C. K. Larive, M. Lin, B. J. Piersma, and W. R. Carper. J.Phys. Chem. , 99:12409–12412, 1995. M. Salanne, C. Simon, and P. A. Madden.