Analysis of the heterogeneous dynamics of imidazolium-based [Tf2N−] ionic liquids using molecular simulation

The complex dynamic behaviour of the imidazolium-based ionic liquids [Cnmim+][Tf2N−], n = 4, 8, 12 is examined at various temperatures and at atmospheric pressure using molecular dynamics simulation. An existing all-atom force field is further optimised in order to attain reasonable agreement with experimental data for transport properties, such as self-diffusivities and viscosities. Dynamical heterogeneity phenomena are quantified through the calculation of the non-Gaussian parameter and the deviation of the self-part of the van Hove correlation function from the expected normal distribution. From this analysis, ions that move faster or slower than expected are detected in the system. These subsets of ‘fast’ and ‘slow’ ions form individual clusters consisting of either mobile or immobile ions. Detailed analysis of the ions’ diffusion reveals preferential motion along the direction of the alkyl tail for the cation and along the vector that connects the two sulphur atoms for the anion. For the longest alkyl tails, the heterogeneity in the dynamics becomes more pronounced and is preserved for several nanoseconds, especially at low temperatures.


Introduction
Over the last decade, ionic liquids (ILs) have been subject to extensive experimental and computational study due to the unique combination of a number of physical properties that renders them ideal candidates for use in a wide range of applications [1], from electrochemistry and catalysis to separation processes and environmental engineering. The great number and chemical diversity of the potential anion and cation combination enables the design of task-specific ILs with controlled properties, as the end-use performance of an IL directly depends on the molecular structure of the constituent ions [2,3].
In this direction, molecular simulation can be a valuable tool for understanding the underlying microscopic mechanisms that govern the macroscopic behaviour of a system and for predicting the properties of interest. The reliability of the calculations from a molecular simulation study largely relies on the ability of the force field to accurately model the intra-molecular and inter-molecular interactions. It is well known that classical force fields with fixed partial charges and a total ionic charge equal to ±1e − underestimate the transport properties of ILs [4][5][6][7]. It has been shown [8] that taking into account polarisability and charge transfer effects yields a more realistic representation of the * Corresponding author. Email: nvergad@chem.demokritos.gr potential energy of ILs and thereby, better predictions of their transport properties. However, many-body polarisable force fields are much more complex to implement and their use significantly increases the computational cost compared to using classical force fields. To address this issue, classical force fields with reduced ionic total charge have been widely used in simulation studies of ILs as an alternative way to incorporate charge transfer and polarisability effects [9][10][11]. A very comprehensive review on this topic is given in Ref. [6].
ILs exhibit very complex dynamic behaviour compared to conventional fluids [12][13][14][15][16][17][18][19][20] and their dynamics resembles that found in polymers, especially at low temperatures where their behaviour is almost glassy like, spanning a wide range of time scales for the relaxation of the various modes in the system. Microscopic heterogeneity phenomena, similar to the ones reported for glass-forming materials and super-cooled liquids [21][22][23] have been observed in ILs both in experiments [23,24] and in molecular simulations. Computational studies on the heterogeneous nature of the dynamics of ILs have focused on the quantification of the non-Gaussian character in the dynamics by calculating the non-Gaussian parameter and by analysing the deviations of the self-part of the van Hove correlation func-C 2014 Taylor & Francis tion [25] from the expected Gaussian form. Del Pópolo and Voth [16] studied the dynamics of [C 2 mim + ][NO 3 − ] and identified spatial correlations between ions of similar mobility. Hu and Margulis [17,18] conducted molecular dynamics (MD) simulations on [C 4 mim + ][PF 6 − ] and investigated the coupling between translational and rotational diffusion and the connection of dynamic heterogeneity phenomena to the 'red-edge effect'. Köddermann et al. [26] examined the connection between dynamic heterogeneity phenomena and the non-validity of the Stokes-Einstein and the Stokes-Einstein-Debye relations by simulating [C 2 mim + ][TF 2 N − ] and its mixtures with chloroform. Liu and Maginn [13] analysed the non-Gaussian behaviour of [C 4 mim + ][TF 2 N − ] at various temperatures and studied the anisotropy in the ions' translational dynamics, while Urahata and Ribeiro [14] focused on the investigation of the link between dynamic heterogeneity and the conformations of alkyl chains in [C 4 mim + ][Cl − ] IL. Mutual diffusion of anions and cations and the multi-fractal nature of the heterogeneous dynamics have been thoroughly studied by Habasaki and Ngai [19].
In this work, the complex dynamics of various 1-methyl-imidazolium bis(trifluoromethylsulfonyl) imide ILs, [C n mim + ][Tf 2 N − ] (n = 4, 8,12), are studied using MD simulation. This is a continuation of our previous study on these systems [20] in which extensive analysis of the effect of the cation alkyl chain length on the thermodynamic, structural, segmental dynamics and transport properties was performed in a wide temperature range and at atmospheric pressure. A wealth of information on the microscopic mechanisms that characterise the spatial organisation and govern the dynamic behaviour of these systems was revealed. Moreover, the need to adequately incorporate polarisability effects in the force field in order to obtain reliable predictions on the transport properties was noted.
In this study, an optimised force field is proposed and new results on densities, ion self-diffusivities and viscosities are reported. The ions' translational motion is analysed into specific axes dictated by their geometry, and anisotropy in the diffusion is detected in both ions. The heterogeneity in the dynamics is thoroughly investigated by calculating the non-Gaussian parameter and the self-part of the van Hove function. Ions that move faster or slower than expected are being detected through deviations from the normal distribution. Dynamically distinguishable ions are found to be strongly spatially correlated. The effect of the temperature and the alkyl chain length on the above properties is examined in depth.

Simulation details
Long MD simulations of 100 [C n mim + ][Tf 2 N − ] (n = 4, 8, 12) ion pairs were performed at 298.15, 348.15 and 398.15 K and atmospheric pressure utilising NAMD software [28] and using as initial configuration for each run the last snapshot of the well equilibrated trajectories (∼80 ns) of our previous study [20]. At first, a 15-25 ns simulation in the isobaric-isothermal ensemble (NPT) was performed in order to relax the system and calculate the density. Subsequently, simulations in the canonical ensemble (NVT) of the order of 20-40 ns -depending on the energy equilibration -were performed at the mean density obtained from the NPT simulations, in order to calculate the transport properties. Snapshots were stored every 1 ps during each simulation. The Lorentz-Berthelot mixing rules were used to calculate the force-field parameters between unlike atoms. A Langevin piston was used to control temperature with a 5 ps −1 damping factor and the Nosé-Hoover barostat was used for the pressure control with oscillation period equal to 200 fs and a damping factor of 100 fs. Electrostatics were handled by means of particle-mesh Ewald method while the reversible reference system propagator algorithm [27,29] was used as a multiple time step algorithm with a 1 fs reference time step in order to speed up the MD simulations. Short-range non-bonded van der Waals interactions were calculated every 2 fs and full electrostatic interactions were computed every 4 fs. A cut-off distance of 12 Å was used to truncate the van der Waals interactions using a switching function beyond 10.5 Å . A long-range correction was applied to the system's energy and virial to account for the neglected van der Waals forces due to switching and cut-off of the Lennard-Jones potential.

Force field
The all-atom force field used to model the potential energy (intra-molecular and inter-molecular interactions) of the system is given by the following expression: where b, θ , χ and ψ denote bond length, bond angle, dihedral and improper angle, respectively, and the subscript '0' refers to the equilibrium values. Parameter n in the term of the dihedral potential is the multiplicity of the dihedral angle while δ is the phase shift of the dihedral potential over the full range of rotation. Partial charges are denoted by q i , ε 0 is the vacuum permittivity and ε, σ are the Lennard-Jones parameters. Electrostatic and Lennard-Jones interactions between atoms separated by less than three covalent bonds were excluded. For atoms separated by three covalent bonds, Coulombic interactions were scaled by a factor of 0.4 while Lennard-Jones interactions were fully included. Non-bonded terms were fully considered for all other interactions. Bonded and Lennard-Jones parameters were obtained from Cadena and Maginn [4] and from Cadena et al. [30] for the cations and Canongia Lopes and Pádua [31] for the anion. The partial charges were derived from quantum mechanical calculations performed on an isolated [C n mim + ][Tf 2 N − ] (n = 4, 8, 12) ion pair using a density functional theory (DFT) methodology. Details on the procedure followed for the calculation of the partial charges can be found in Refs [20] and [32]. The calculation of the charges based on the minimum energy conformations of the ionic pair, results in reduced total ions' charge, a fact that is attributed to the charge transfer and polarisation effects [6,7]. The results on thermodynamic, structure and transport properties of the three ILs using directly the charges from the quantum mechanical calculations are reported in Ref. [20]. The thermodynamic properties were in good agreement with the experimental data, whereas there was a systematic underestimation and overestimation by an order of magnitude of the self-diffusion coefficients and the viscosity, respectively. In the present study, the charge distribution was uniformly scaled to a total ionic charge equal to ±0.75e − for all three ILs. Reduced total ionic charges have been widely used in molecular simulation studies of ILs [6,9,13,33] and have been reported to result in better agreement with experimental data on transport properties. In our analysis, the uniform scaling of the partial charges results in significant improvement of the predicted transport properties, whereas the density is very weakly affected and there is no notable impact in the structural properties. Calculations performed with the force field used in our previous work will be referred to as '6RC' since the charge distribution was derived from ab initio calculations on six minimum energy conformations of the ion pair, while results using the new force field will be referred to as '6RC scaled 0.75e'. All the force-field parameters used in this work are provided as Supplementary Information (Table S1-S6).

Density, self-diffusivities and viscosity calculations
MD predictions for the density of the three ILs are plotted in Figure 1 against experimental data [3,[34][35][36] as a function of temperature at 1 bar. Our previous calculations are also shown for comparison (6RC). In all cases, the deviation between the calculated values and the experimental measurements is less than 1%. There is a systematic decrease of the values compared to our previous calculations, primarily driven by the reduction in the intensity of the at-  Table S7 in the Supplementary Information. For the calculation of the ions' self-diffusion coefficient, the Einstein relationship was used, where r i (t) is the position of the ith ion's centre of mass at time t and the brackets denote the mean-square displacement (MSD) over all ions' centres of mass. The diffusion coefficients were calculated in the Fickian regime, which is identified by a slope equal to unity in the log(MSD) versus log(t) plot.  Figure 2(a)-(c), respectively, as a function of temperature, together with experimental data [35,37]. For comparison, our previous results on self-diffusivity [20] are also plotted (6RC). There is a clear improvement in the accuracy of the calculations with the new force field. Simulation captures quite well the temperature dependence of the self-diffusivities as well as the fact that the cation's self-diffusion coefficient is higher than the anion's for [C 4 mim + ][Tf 2 N − ]. This tendency is also clear in the available experimental data that are plotted in Figure 2(a) and has been reported in past simulation studies of imidazolium-based ILs [7,11,13,15,38,39]. As the alkyl tail becomes longer, the cation's mobility is reduced and the ions attain comparable self-diffusivities, as can be seen in Figure 2 three ILs are provided in Table S8 in the Supplementary Information. Calculations of the shear viscosity were performed using the Green-Kubo relationship: where P ij (t) is the ij-element of the pressure tensor at time t (i =j). Shear viscosity is a collective property of the system and, therefore, its calculation is prone to errors due to the sensitivity of the stress correlation function to the intramolecular vibrational modes. At the same time, experimental measurements of viscosity are also very sensitive to impurities of the testing material; consequently, even experimental results obtained using the same method deviate from each other [40]. The stress correlation function was calculated at the last section of each productive NVT run during which the information of the pressure was stored every 8 fs and the viscosity was estimated by the integral plateau value after reassuring that the pressure autocorrelation function was fluctuating around zero. In Figure 3 [20] and are compared to experimental data [35,[41][42][43]. The simulations using the optimised force-field result in a much better agreement between the calculated viscosity values and the experimental data, and the temperature dependence is captured very accurately. The calculated viscosity values are provided in Table S9 in the Supplementary Information.

Analysis of the ions' translational motion
In an effort to investigate the origin of the fact that in many cases cations are reported to have enhanced diffusivity compared to anions, Urahata and Ribeiro [15]  studied the anisotropy in the cation's diffusivity for the case of imidazolium-based ILs with small-sized anions. They analysed the translational motion of the cation along specific axes, namely along the vector that is normal to the imidazolium ring, along the vector NN that connects the two nitrogens of the imidazolium ring and along the vector NNp that is perpendicular to the other two vectors. They identified an enhanced movement of the cation along the NNp vector that lies on the plane of the imidazolium ring and is almost perpendicular to the alkyl chain. Liu and Maginn [13] performed a similar analysis for the case of [C 4 mim + ][Tf 2 N -] and also found anisotropy in the diffusional motion, detecting in their case the direction along the NN vector as the most favourable one.
We further investigated this phenomenon and examined the effect of the alkyl tail length by analysing the cations' translational motion along the three axes mentioned above and, additionally, along the N3-Ct vector that connects the terminal carbon atom of the alkyl tail (Ct) and the N3 nitrogen of the imidazolium ring (see Figure 4(a)). At the same time, the anions' translational motion was analysed along the axes: (1) SS that connects the two sulphurs, (2) SSp that is perpendicular to SS and lies on the plane defined by the two sulphur atoms and the nitrogen of the anion and (3) the vector normal to the plane that is defined by the SS and the SSp vectors (see Figure 4(b)). The MSD of the ions' centre of mass in each specific axis is compared to the 1/3 of the total MSD of the ions' centre of mass.
In Figure 5, the MSD along each axis is shown for the cation and the anion of [C 4 Figure 5(c)), the MSD along N3-Ct at 8 ns is almost 50% higher than the 1/3 of the total MSD, whereas for the case of [C 12 mim + ][Tf 2 N − ] that has the longest alkyl tail, this deviation reaches almost 100% ( Figure 5(e)). This behaviour is preserved at very long times and is present even at time scales that correspond to the Fickian regime of the total MSD at 298.15 K. The analysis of the anion's diffusion also reveals a favourable motion along the SS vector, which is more pronounced for the two ILs with the longest alkyl tails (see Figure 5(b), 5(d) and 5(f)). As the temperature increases, the difference in the MSDs is only present at very short time scales, for example, for [C 4 mim + ][Tf 2 N − ] at 398.15 K, this preferential translational motion is observed below 1 ns for the cation and below 0.5 ns for the anion (see Figure S1 in the Supplementary Information).

Dynamic heterogeneity
Ions' translational motion is further analysed in order to investigate the presence of heterogeneities in the dynamics by detecting deviations from the expected Gaussian behaviour. The global non-Gaussian behaviour of a system can be quantified by analysing the time dependence of the non-Gaussian parameter a 2 (t) defined as [44] where |r i (t) − r i (t 0 )| is the displacement of an ion's centre of mass at time t with respect to its position at time t 0 . In the case of a Gaussian distribution of |r i (t) − r i (t 0 )|, a 2 (t) is equal to zero, while any deviation from zero is indicative of dynamic heterogeneities in the system revealing simultane-ously the time scale at which the non-Gaussian behaviour is present. The time at which a 2 (t) reaches a maximum value indicates the time interval of the most intense dynamic heterogeneity.
The self-part G s (r, t) of the van Hove correlation function [25] is given by and corresponds to the probability that a particle (or the centre of mass of a molecule) is at position r i (t) at time t given that it was located at r i (0) at time 0. Note that for all times t: The probability density that a particle moves by distance r is given by the solid angle integral [45]: In case of isotropy, G s (r, t) = G s (r, t) and g s (r, t) = 4πr 2 G s (r, t).
When heterogeneous dynamics is present, the self-part G s (r, t) of the van Hove correlation function deviates from the expected Gaussian distribution that is given by the following expression: where | r(t)| 2 is the MSD of the particles. In Figure 6, the non-Gaussian parameter a 2 (t) is plotted against time for the cation and the anion of the three ILs at 298.15, 348.15 and 398.15 K. For both ions, a 2 (t) maximises at intermediate times and as the temperature decreases, the maximum of a 2 (t) tends to appear at longer times reaching a higher value. Cations systematically exhibit a higher degree of heterogeneity than anions that is reflected to a higher peak of a 2 (t). These findings are in agreement with previous analysis of dynamical heterogeneities of ILs in the literature [13,16]. It is evident, especially at low temperatures, that as the alkyl tail becomes longer, the heterogeneity in the dynamics is present for much longer time intervals for both ions since the non-Gaussian parameter starts to decay at much longer times. This indicates the existence of a cage regime of longer duration, a fact that can be directly related to the MSD behaviour and the time scales involved for reaching in each case the Einstein regime of diffusion. For example, in the case of the behaviour of the non-Gaussian parameter of the cation of [C 12 mim + ][TF 2 N − ] at 298.15 K (Figure 6(e)), the a 2 (t) starts to decay at about 8 ns, that is, the time at which the cation's diffusion reaches the Fickian regime (see Figure S2 in the Supplementary Information).
The self-part of the van Hove function 4πr 2 G s (r, t) calculated at t = 10, 300, 1000 and 3000 ps at 298.15 K is plotted in Figure 7 as a function of distance r and in comparison to the corresponding expected Gaussian distribution These plots clearly show that there are deviations from the Gaussian distribution and that heterogeneity phenomena are present at time intervals from a few ps up to several ns depending on the IL. This is more pronounced for the two ILs with the longest alkyl tails, for which the caging regime involves longer times. In agreement with the a 2 (t) results, the G s (r,t) of the cations exhibit greater deviations from the Gaussian distribution, a fact that is attributed to the non-uniform way that anions surround the cations, while the anions experience a more uniform neighbourhood [13]. In all cases, the deviations from the Gaussian distribution are reduced as the temperature increases (see for example, Figure S3 in Supplementary  Information).
Following the approach of Kob et al. [21] on supercooled liquids, ions can be classified according to their mobility (the distance each ion has travelled at time t). In many cases G s (r, t) and G g (r, t) intersect both at short, r 1 , and at long distances, r 2 , indicating the existence of ions with low and high mobility, i.e. G s (r, t) > G g (r, t) for r < r 1 and G s (r, t) > G g (r, t) for r > r 2 , respectively (see for example Figure 7(e)). The spatial correlation between the groups of ions with different mobility is examined by calculating all the relevant pair distribution between anion and cation centres of mass in the bulk. The term 'bulk' refers to the total radial distribution function calculated without any discrimination on the ions' mobility. The subsets of ions considered in the calculations of g(r) consist of the 8%-10% most immobile (if present) and mobile ions. These plots clearly depict that there is strong spatial organisation between anions and cations of the same mobility (fast-fast or slow-slow) as indicated by the higher peak of the g(r) in comparison with the g(r) in the bulk, while   (Figure 8(a)) and as the alkyl tail becomes longer the peak becomes wider but in all cases it is observed at the same distance with the corresponding ones of the bulk. Similar behaviour is observed for slow anion and slow cation pair correlation functions for the [C 8 mim + ][Tf 2 N − ] (Figure 8 thermore, the anion-anion pair distribution function shows similar characteristics: fast-fast and slow-slow anion correlation is more intense that the one in the bulk, while it is less likely to find anions of different mobility close to each other, as shown in Figure 9. These observations are much weaker when calculating the corresponding cationcation correlation functions as depicted in Figure 10. The time evolution of the spatial organisation of the mobile cations and anions is shown in Figure 11

Conclusions
An existing all-atom force field for [C 4 mim + ][Tf 2 N − ], [C 8 mim + ][Tf 2 N − ] and [C 12 mim + ][Tf 2 N − ] ILs was further optimised in order to attain good agreement with the experimental data on transport properties by uniformly scaling the charge distribution to yield a ±0.75e − total ion charge. Long MD simulations of several tens of nanoseconds were performed in order to calculate the transport properties of the three ILs mentioned above and analyse their complex dynamic behaviour focusing on the investigation of microscopic heterogeneity phenomena. The results of the optimised force field on density, ion self-diffusion coefficients and viscosity exhibit a significant improvement in the agreement with the experimental data. Detailed analysis of the ions' translational motion reveals anisotropy in the ions' diffusion: there is a preferential movement along the direction of the alkyl tail for the cation and along the direction of the vector that connects the two sulphurs for the anion. The heterogeneous nature of the dynamics is quantified by the non-Gaussian parameter and is further analysed by comparing the self-part of the van Hove function to the expected Gaussian distribution. For the longest alkyl tails, the heterogeneity in the dynamics becomes more pronounced and is preserved for several nanoseconds at low temperatures. Subsets of 'fast' and 'slow' ions are defined by the deviation between the van Hove function and the Gaussian distribution. Dynamically distinguishable ions are found to be strongly correlated in space forming clusters of ions of the same mobility.

Funding
The MD simulations have been largely performed under the HPC-EUROPA2 project (project number 228398) with the support of the European Commission -Capacities Area -Research Infrastructures. Financial support from the Seventh European Framework Programme for Research and Technological Development for the project 'Novel Ionic Liquid and Supported Ionic Liquid Solvents for Reversible Capture of CO 2 ' (IOLICAP Project 283077) is gratefully acknowledged. This work is also part of the Project "THALIS" which is implemented under the Operational Project "Education and Life Long Learning" and is co-funded by the European Union (European Social Fund) and National Resources (ESPA).

Supplemental data
Supplemental data for this article can be accessed here.