Molecular dynamics simulation of geminal dicationic ionic liquids [C_n(mim)₂][NTf₂]₂ – structural and dynamical properties

Abstract

In this work, the structural and dynamical properties of two imidazolium-based geminal dicationic ionic liquids (GDILs), i.e. [C_n(mim)₂][NTf₂]₂ with n = 3 and 5, have been studied to obtain a fundamental understanding of the molecular basis of the macroscopic and microscopic properties of the bulk liquid phase. To achieve this purpose, molecular dynamics (MD) simulation, density functional theory (DFT) and atoms in molecule (AIM) methods were used. Interaction energies, charge transfers and hydrogen bonds between the cation and anions of each studied GDIL were investigated by DFT calculations and also AIM. The mean square displacement (MSD), self-diffusion coefficient, and transference number of the cation and anions, and also the density, viscosity and electrical conductivity of the studied GDILs, were computed at 333.15 K and at 1 atm. The simulated values were in good agreement with the experimental data. The effect of linkage alkyl chain length on the thermodynamic, transport and structural properties of these GDILs has been investigated. The structural features of these GDILs were characterized by calculating the partial site–site radial distribution functions (RDFs) and spatial distribution functions (SDFs). The heterogeneity order parameter (HOP) has been used to describe the spatial structures of these GDILs and the distribution of the angles formed between two cation heads and the middle carbon atom of the linkage alkyl chain was analyzed in these ILs. To investigate the temporal heterogeneity of the studied GDILs, the deviation of the self-part of the van Hove correlation function, G_s(r⃑,t), from the Gaussian distribution of particle displacement and also the second-order non-Gaussian parameter, α₂(t), were used. Since, the transport and interfacial properties and ionic characteristics of these GDILs were studied experimentally in our previous studies as a function of linkage chain length and temperature, in this work, we try to give a better perspective of the structure and dynamics of these systems at a molecular level.

---

1. Introduction

Ionic liquids (ILs) are molten salts having melting points less than 100 °C.¹ ILs show specific properties such as negligible vapor pressure, low melting points, non-flammability, high thermal stability and high intrinsic ionic conductivity which make them suitable for industrial applications.² Due to their unique properties, ILs are described as potential “green” solvents to replace conventional organic ones to reduce environmental pollution.³ The possibility of selection of different cations and anions in the structure of ILs allows for a large variety of tunable properties and applications.⁴

Geminal dicationic ionic liquids (GDILs), as a new category of ILs, consist of two head groups linked by a rigid or flexible spacer and two monoanions.⁵ These compounds have received considerable attention from researchers in the academic and industrial fields due to their higher melting points, wider liquid range, and better thermal stability compared with the conventional monocationic ionic liquids (MILs). They can be used as solvents in organic reactions,⁶ as lubricants,^7–10 as the stationary phase in gas chromatography columns,^11–13 and as electrolytes in secondary batteries,¹⁴ in electrospray ionization mass spectrometry¹⁵ and in dye sensitized solar cells especially at high temperatures.^16,17 Although, the number of reports on DILs is increasing gradually, the studies on this group of ILs are still limited and the knowledge about them is much less than that of MILs. There are few experimental studies on these compounds to understand their structural and dynamical properties and also to explain these properties at a molecular level.^8,18–21 For example, Majhi and Sarkar²² investigated the microscopic structural organization of neat ILs and IL-based gels experimentally through resonance energy transfer (RET) studies. A small number of computational studies on imidazolium-based GDILs have been reported so far.^4,21,23–27 Bodo et al.²⁷ studied long range structural ordering of [C_n(mim)₂] [NTf₂]₂ (n = 3, 6, 9 and 12) and compared their bulk liquid structures with gas-phase structures and also with the structures of the monoimidazolic counterparts. Yeganegi et al.²⁶ investigated the effect of the anion type and the alkyl linkage length on the diffusion, density and microscopic structure of the liquid phase for (C_n(mim)₂X₂) (n = 3, 6, or 9) where X = PF₆⁻, BF₄⁻, and Br⁻. Li et al.²⁵ studied the effects of the linkage chain length on the nanoscale organization of [C_n(mim)₂](BF₄)₂ (n = 3, 6, 9, 12, 16) compared with those of the free alkyl chain length in MILs. They proposed idealized structural models for nanoaggregates formed by the alkyl chains of DILs and MILs. Also, Lovelock et al.²⁴ studied the vaporization of [C₃(mim)₂][NTf₂]₂ by means of the ioncyclotron resonance mass spectroscopy technique and MD simulation. However, more studies on some features of these compounds, such as microscopic interionic interactions, heterogeneity, and their spatial structures at a nanoscopic level seem to be of great importance.

In our previous studies,^28,29 we studied experimentally the transport and interfacial properties and ionic characteristics of imidazolium-based short alkyl chain length DILs ([C_n(mim)₂][NTf₂]₂) as a function of linkage chain length and temperature. In the present work, to obtain a fundamental understanding of the molecular basis of the macroscopic and microscopic properties of the bulk liquid phase in the studied GDILs, i.e. [C_n(mim)₂][NTf₂]₂ with n = 3 and 5, in which the cations are 1,3-bis(3-methylimidazolium-1-yl)alkane (propane or pentane) ([C_n(mim)₂]²⁺), and the anion is bis(trifluoromethylsulfonyl)imide ([NTf₂]⁻), we have used molecular dynamics (MD) simulation and also density functional theory (DFT) and atoms in molecule (AIM) methods. Fig. 1 shows the schematic structures and the atomic labels of the cations and anions of the studied GDILs. The effect of linkage alkyl chain length on the thermodynamic, transport and structural properties of these GDILs has been investigated. In addition, the heterogeneity order parameter (HOP) has been used to describe the spatial structures of these GDILs and the distribution of the angles formed between two cation heads and the middle carbon atom of the linkage alkyl chain was analyzed in these ILs. To investigate the temporal heterogeneity of the studied GDILs, we try to study the deviation of the self-part of the van Hove correlation function, G_s(r⃑,t), from the Gaussian distribution of particle displacement and also the second-order non-Gaussian parameter, α₂(t).

Fig. 1 Schematic structures and atomic labels of cations and anions of the studied GDILs. The hydrogens connected to the carbons of the alkyl chains and also CS have not been shown in this figure.

2. Computational details

The initial configuration of each IL (a single ion triplet) has been created as follows: the structures of the isolated cation and anions were optimized by density functional theory (DFT) at the B3LYP/6-31++G(d,p) level of theory, using the GAUSSIAN 03 program.³⁰ The structure of each ion was checked by harmonic vibrational frequencies to ensure the absence of negative frequencies and verify the existence of a true minimum. The charges of the separate cation and anions were set at +2e and −e, respectively. After that, the structure of each IL (cation and anions together) was optimized using the B3LYP/6-31G(d) level of theory. All atom non-polarizable force fields were employed to simulate the studied DILs, which were taken from the work of Tsuzuki et al.³¹ The potential energy of the systems was modeled using the following standard functional form:

The functional form of the force field contains four kinds of potential energies: stretching of covalent bonds, bending of angles, torsion around dihedral angles, and nonbonded interactions. In the first three terms, the intermolecular bonds and angles are described by harmonic terms and dihedral torsion energies are represented by a series of cosines. Nonbonded interactions are described by the last term of eqn (1), including the Lennard-Jones (12-6) potential for the van der Waals interactions and the Coulombic term for the electrostatic interactions between partial charges placed on the atomic sites. The cross term parameters i.e. ε_ij and σ_ij are obtained by combining and mixing rules, i.e. ε_ij = (ε_iiε_jj)^1/2 and σ_ij = (σ_ii + σ_jj)/2, respectively. The atomic partial charges were calculated using the CHarges from ELectrostatic Potentials using a Grid (CHELPG) method³² at the B3LYP/6-31++G(d,p) level of theory. The total charges for the ions obtained in this work are considerably less than unity which are in agreement with the previous simulation studies on other IL systems.^33,34 The electrostatic charges, bond lengths, bond angles and dihedrals for the studied GDILs were taken from our ab initio calculations. The electrostatic charges have been reported in the ESI† (Tables S1 and S2).

The binding energy for each ionic triplet (E_b) was defined as:

E_b = E^gCP-D3_GDIL − (E^gCP-D3_cation + 2E^gCP-D3_anion) (2)

where E^gCP-D3_cation, E^gCP-D3_anion, and E^gCP-D3_GDIL are the corrected total energy of the cation, anion, and GDIL, respectively. The values of E_b were calculated as the difference among the total energy of each GDIL and the total energies of its constituent ions optimized at the same level of theory (B3LYP/6-31G(d)). The total electronic energies for each cation and anion, and for the ionic triplet, were corrected using the Grimme's DFT-gCP-D3 method.^35,36E^gCP-D3 for each configuration was obtained using eqn (3).

E^gCP-D3 = E^el + E^gCP + E^D3 (3)

where E^el is the electronic energy and E^gCP and E^D3 stand for corrections of geometrical counterpoise energy and dispersion, respectively. To get a deeper insight into the anion/cation interactions and bond characteristics, quantum theory of atoms in molecule (QTAIM) was used via the AIM2000 program package.^37,38

All of the ILs were simulated using the DL_POLY 2.17 molecular dynamics package.³⁹ For each IL, a periodic cubic box containing 216 ion triplets was used. In these boxes, the total number of atoms of [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂, was 13 608 and 14 904, respectively. Temperature and pressure were controlled using the Nosé–Hoover thermostat and barostat⁴⁰ with 0.1 and 2 ps relaxation times, respectively. Long-range electrostatic interactions were calculated with the Ewald summation technique with a precision of 1 × 10⁻⁵. A cut-off distance of 25.0 Å was chosen for the Lennard-Jones interactions and cubic periodic boundary conditions were applied. The equations of motion were solved using the Verlet–Leapfrog integration algorithm.⁴¹ Due to the slow dynamics of these systems, in each case, the simulation was started at a relatively high temperature (350 K) at ambient pressure (1 atm). Equilibration was carried out following an annealing schedule. The temperature of the systems was increased at intervals of 50 K up to 600 K. At each temperature, the simulation was performed for 500 ps under constant NPT conditions at 1 atm. The system was then cooled down sequentially to the target temperature at steps of 50 K. After equilibration, the lengths of the sides of the simulation boxes were 55.10 Å and 57.57 Å for [C₃(mim)₂][NTf₂]₂, and [C₅(mim)₂][NTf₂]₂, respectively. Then, a 2 ns NPT simulation was carried out in order to establish the density and structure of each IL at a specified temperature at ambient pressure. The transport properties of the ions of each GDIL were calculated at 333.15 K from a 8 ns production run under NPT conditions.

The reliability of molecular dynamics simulations depends on the ability of the parameterized force fields to reproduce accurately relevant thermophysical properties such as liquid density. The simulated densities for the studied GDILs are in good agreement with the experimental data and deviations between simulated and experimental densities are about 3% over the whole temperature range (see Table S3, ESI†).

3. Results and discussion

3.1. Quantum chemistry calculations

Calculating the charge transfer upon ionic triplet formation is an important factor for the analysis of the studied GDILs. While polarizable force fields are more successful in predicting the dynamical properties of ILs, they need a noticeably higher computational cost compared with the nonpolarizable ones.⁴² Although, the nonpolarizable force fields greatly underestimate the ionic dynamics, they can reproduce the thermodynamic and structural properties of ILs with a good accuracy.⁴³ On the other hand, when all-atom force fields are used, the computed dynamical properties of ILs converge slowly and the use of polarizable force fields for simulation of ILs is not feasible. To improve the nonpolarizable force fields, we can calculate charge transfers between ions and even use the scaling factor to reduce the partial electrostatic charges on all interaction sites. It has been shown that the use of a scaling factor is an adequate approximation to reproduce the experimental data of ILs, especially for their transport properties within the framework of classical MD studies.^43,44 As mentioned in Section 2, in this work, charge transfer processes between ions were studied through the CHELPG method. This method showed that the net charges on the cations are +1.55e and +1.65e for [C₃(mim)₂][NTf₂]₂, and [C₅(mim)₂][NTf₂]₂, respectively (which means that the net charge on each anion would be approximately −0.77e and −0.82e, respectively). As it can be seen, there is remarkable charge transfer between the corresponding ions which is in agreement with the literature for other types of ILs.^33,45,46

The intermolecular interactions were characterized through the estimation of binding energies as well as AIM analysis. The total electronic energies (E_el), gCP correction energies (E_gCP), D3 correction energies (E_D3) and binding energies are tabulated in Table 1. The reported results show hydrogen bonding between HA hydrogen atoms of the rings with nitrogen or oxygen atoms of the [NTf₂]⁻ anions and also hydrogen atoms of the methylene groups connected to the rings with oxygen of the anions. In agreement with the shorter hydrogen bond distances for most of the presented hydrogen bonds in [C₃(mim)₂][NTf₂]₂ compared with [C₅(mim)₂][NTf₂]₂, the interaction between the [C₃(mim)₂]²⁺ cation and [NTf₂]⁻ anions is stronger than that for [C₅(mim)₂]²⁺ with the anions (Fig. 2). The values of interaction energies between the dication and two anions in the studied GDILs were computed to be −927.402 kJ mol⁻¹ and −832.842 kJ mol⁻¹ for [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂, respectively. These values are much greater than the reported interaction energies between ion pairs in MILs.^4,47,48 This can be due to the stronger electrostatic attractions between the cation and the anions in GDILs compared with the interactions between the ions in MILs which is responsible for the higher melting point and thermal stability of GDILs. This result is in agreement with those obtained by Sun et al. in which the interaction energy in [C₃(mim)₂][Br]₂ was more than double the interaction energy in [C₄mim][Br].⁴

Table 1 Electronic energy, gCP and D3 corrections and binding energies of GDILs and their corresponding ions

Configuration	E ^el (a.u.)	gCP correction (a.u.)	D3 correction (a.u.)	E ^el-gCP-D3 (a.u.)	E b (kJ mol^−1)
[C3(mim)2]^2+	−648.498	0.080	−0.055	−648.473	—
[C5(mim)2]^2+	−727.150	0.095	−0.064	−727.119	—
[NTf2]^−	−1827.205	0.103	−0.034	−1827.136	—
[C3(mim)2][NTf2]2	−4303.267	0.326	−0.158	−4303.099	−927.402
[C5(mim)2][NTf2]2	−4381.875	0.329	−0.164	−4381.709	−832.842

---

Fig. 2 Structures for [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂ ionic triplets obtained from gas phase calculations at the B3LYP/6-31G(d) theoretical level. ΔE stands for counterpoise corrected interaction energies. Dashed lines show possible hydrogen bonds. Atom color code: gray, carbon; blue, nitrogen; red, oxygen; and white, hydrogen.

To gain a better insight into the intermolecular interactions between the cation and anions in the studied GDILs, AIM analysis was performed on the optimized geometries. In AIM, there is one bond critical point (BCP) between each pair of atoms bonded or interacting, where a critical point is defined by the point localized between two attractors. The AIM theory, based on a topological analysis of electron density at the bond critical point (ρ_BCP) and its laplacian (∇²ρ_BCP), provides a universally applicable tool for the classification of the bonding interactions occurring in any molecular system. The ρ_BCP is used to describe the strength of a bond, with stronger bonds associated with a larger ρ_BCP value. The ∇²ρ_BCP describes the characteristic of the bond. When ∇²ρ_BCP < 0, it is named as a covalent bond. When ∇²ρ_BCP > 0, it refers to a closed-shell interaction⁴⁹ with characteristics of ionic bonds, H-bonds, or vdW interactions. In this paper, the values of ρ_BCP and ∇²ρ_BCP have been observed only for intermolecular H-bonds. The considered criterion to define a hydrogen bond according to the AIM approach is that ρ_BCP and ∇²ρ_BCP must be in the 0.002–0.035 a.u. and 0.024–0.139 a.u. ranges, respectively.⁵⁰ For both of the studied GDILs, the values of ρ_BCP and ∇²ρ_BCP which are given in Table 2 are in the above mentioned ranges. Also, the topological graphics of the lowest energy ionic triplets of the studied GDILs are shown in Fig. 3. In this figure, several interaction sites can be determined between the cation and anions in the studied GDILs, where the anion interaction sites are nitrogen and oxygen atoms and the cation interaction sites are HA hydrogen atoms of the rings and also hydrogen atoms of the methylene groups connected to the rings. In this figure, the most important possible hydrogen bonds are shown with green lines. The values of ∇²ρ_BCP listed in Table 2 are all positive, indicating the existence of the electrostatic characteristics in the inter-ionic hydrogen bonds. According to Table 2, in both of the studied GDILs, the interaction of oxygen atoms of the anions with the HA atom of each ring of a cation is stronger than the interaction of nitrogen atoms of the anions with HA. Furthermore, it can be seen that the interaction of oxygen atoms of the anions with an HA atom of each ring is much greater than the interactions with the hydrogens of methylene groups connected to the rings. Thus, the charge depletions between the HA atoms and also hydrogens of methylene groups from the cation and the oxygen and nitrogen atoms of the anions correspond to strong hydrogen bonds.

Table 2 Electron density (ρ) and Laplacian of the electron density (∇²ρ) in the BCPs (bond critical points) for the studied GDILs

Structure	Bond	ρ BCP	∇^2ρBCP
[C3(mim)2][NTf2]2	N(anion)⋯HA	0.0162	0.0519
	N(anion)⋯HA	0.0161	0.0518
	O(anion)⋯HA	0.0164	0.0548
	O(anion)⋯HA	0.0165	0.0550
	O(anion)⋯H(Cs)	0.0128	0.0459
	O(anion)⋯H(Cs)	0.0128	0.0458
[C5(mim)2][NTf2]2	N(anion)⋯HA	0.0156	0.0548
	N(anion)⋯HA	0.0172	0.0563
	O(anion)⋯HA	0.0195	0.0593
	O(anion)⋯HA	0.0170	0.0557
	O(anion)⋯H(Cs)	0.0097	0.0337
	O(anion)⋯H(Cs)	0.0090	0.0318

---

Fig. 3 Topological graphics of the lowest energy ionic triplets obtained from AIM analysis of the (a) [C₃(mim)₂][NTf₂]₂ and (b) [C₅(mim)₂][NTf₂]₂ showing bond critical points (BCP, red circles), ring critical points (yellow circles), and cage critical points (green circles). For the sake of visibility, only bond paths (pink lines) have been shown. The most important possible hydrogen bonds have been shown with green lines.

3.2. Molecular dynamics simulations

3.2.1. Structural properties. To obtain a better insight into the local structure and organization of the studied GDILs, the site–site and center of mass radial distribution functions (RDFs) which give the probability of finding a pair of atoms at a given distance were computed. Fig. 4 shows the analysis of the site–site RDFs for the oxygen and nitrogen atoms of the anion with the hydrogen atoms of the cation rings, i.e. HA, H₁ and H₂ for both of the studied GDILs at 333.15 K and 1 atm. The peaks due to the interactions between oxygen atoms of the anions and the hydrogens of the cation rings exhibit a strong first peak at about 3 Å, indicating a preferential interaction over the interactions between the nitrogen atom of the anion and the hydrogen atoms of the rings, which exhibit first peaks at about 5 Å. In other words, the first peaks of the interactions between a nitrogen atom and the hydrogens of a ring are positioned at larger separations. The peaks of oxygen atoms of the anions with all of the ring hydrogen atoms are more intense and appeared at shorter distances. These RDFs show that, compared with the nitrogen atoms, the oxygen atoms of the anions are more likely to form hydrogen bonds with the ring hydrogens. This result is in agreement with the shorter length of hydrogen bonds of oxygens with the ring hydrogens in DFT and AIM calculations. In each case of interactions between oxygen or nitrogen atoms of the anions with the hydrogen atoms of the cation ring, the peak intensity of HA is stronger than those of H₁ and H₂. Fig. 4 also shows that the interactions between the oxygen and nitrogen atoms of the anion with hydrogen H₂, which is toward the inside of the cation, is slightly stronger than the interactions with H₁, which is toward the outside of the cation.

Fig. 4 Radial distribution functions of the oxygen and nitrogen atoms of the anion with the ring hydrogen atoms, i.e. HA, H₁ and H₂ for (a) [C₃(mim)₂][NTf₂]₂ and (b) [C₅(mim)₂][NTf₂]₂ at 333.15 K and 1 atm.

The radial distribution functions of the oxygen atoms of the anions with the hydrogen atoms of the methyl groups (hydrogen atoms connected to the CS atom) for the studied GDILs are given in Fig. 5 at 333.15 K and 1 atm. As this figure shows, the first peaks of the RDFs show less intensity than those observed for the interactions of oxygen atoms with the cation ring hydrogens, i.e. HA, H₁ and H₂ (see Fig. 4). This interaction is a little stronger for [C₃(mim)₂][NTf₂]₂, which is in agreement with higher values of ρ_BCP obtained from AIM calculations.

Fig. 5 Radial distribution functions of the oxygen atoms of the anion with the hydrogen atoms of the methyl groups (hydrogen atoms connected to CS atom) for (a) [C₃(mim)₂][NTf₂]₂ and (b) [C₅(mim)₂][NTf₂]₂ at 333.15 K and 1 atm.

Fig. S1 (ESI†) shows the center of mass RDFs of cation–cation, anion–anion, and cation–anion for the studied GDILs at 333.15 K and 1 atm. While the anion–anion and cation–cation interactions show first peaks at about 8 or 9 Å, the cation–anion RDF exhibits a strong first peak at about 4 Å for each of the studied GDILs. Positioning of these interactions at shorter separations indicates preferential cation–anion interactions over the two other ones. In fact, to preserve the conditions of local electroneutrality, the ions of a given type prefer to be surrounded by ions of opposite charge. All RDFs exhibit a second less-intense but broader peak and also a weak third peak at longer distances, which is due to the long-range Coulombic interactions. The cation–cation RDFs exhibit a broad first peak, with a low intensity compared to those of the other two pair types.

The radial distribution functions of the anions [NTf₂]⁻ around the CR carbon atom of the rings and also the carbon atom in the middle of the alkyl chain for the studied GDILs have been given in Fig. S2 (ESI†) at 333.15 K and 1 atm. The labels are defined according to those in Fig. 1. The distinct and sharp peaks at about 6 Å indicate strong interactions between the mentioned carbon atoms of the cations and the anions. Fig. S2 (ESI†) shows that the interactions between the anions and the CR carbon of the rings are stronger than the interactions of the anions with the carbon atom in the middle of the alkyl chain for the studied GDILs.

Alkyl chain–alkyl chain RDFs of the studied GDILs are shown in Fig. 6 at 333.15 K and 1 atm. As this figure shows, in spite of [C₃(mim)₂][NTf₂]₂, which has a weak broad peak at about 10 Å, [C₅(mim)₂][NTf₂]₂ shows a sharp and intense peak at about 5 Å. This shows the aggregation of the alkyl chains in [C₅(mim)₂][NTf₂]₂. Two kinds of interactions i.e. Coulombic interactions between the rings and anions and the short-range attractive van der Waals interactions between the alkyl chains compete with each other under a geometrical constraint between the rings and alkyl chain groups in these GDILs. It seems that, with increasing the linkage alkyl chain length in [C₅(mim)₂][NTf₂]₂, the rings have less effect on the aggregation of the alkyl chain groups and the alkyl chain groups have more freedom to be closer to each other. In the case of [C₃(mim)₂][NTf₂]₂, the geometrical constraint restricts the aggregation of the alkyl chains due to the short linkages.

Fig. 6 Radial distribution functions of alkyl chain–alkyl chain of the studied GDILs at 333.15 K and 1 atm.

To give useful visual insights into the spatial arrangement between the anions and cations in the studied GDILs, the three-dimensional spatial distribution functions (SDFs) have been calculated via the Travis package⁵¹ and visualized by the VMD program⁵² as shown in Fig. 7. In this figure, mauve stands for [NTf₂]⁻, cyan for carbon, blue for nitrogen, and white for hydrogen. As it can be seen in Fig. 7, the probability of finding anions around the imidazolium rings is higher than that around the alkyl chains and the anions are clustered around the rings. For better presentation and clarity, the clustering of anions around the rings has been presented just for one of the rings in Fig. 7. This result is similar to the results obtained for MILs³³ and some other DILs.²⁶ As it is recognizable in the SDFs, the GDILs bend to maximize the electrostatic interactions of the imidazolium rings with the anions and also to minimize the interactions between ions with the same charges.

Fig. 7 Spatial distribution functions (SDFs) of the anion [NTf₂]⁻ around the cations: Mauve, [NTf₂]⁻; cyan, carbon; blue, nitrogen; white, hydrogen.

To investigate the structural nano-organization of the studied GDILs, the distribution of the angles formed by two vectors pointing from the rings to the center atom of the linkage chain for the studied GDILs is shown in Fig. 8 at 333.15 K. Li et al.²⁵ showed that in spite of MILs in which a micelle-like nanoaggregate model proposed by Triolo et al.^53,54 is acceptable, in DILs, nanoscale organization is more complicated due to the constraint of the linkage alkyl chain between two imidazolium rings and reduction of the flexibility. They showed that in a series of [C_n(mim)₂][BF₄]₂ with different values of n, neither the straight chain model, nor the folded chain model exist and the nanoscale aggregates consist of a mixture of a majority of unfolded and a small portion of highly folded DILs. Fig. 10 shows that there is a relatively wide distribution of angles formed between the rings and the center atom of the linkage chain for both of the studied GDILs. While in [C₅(mim)₂][NTf₂]₂, the most probable distribution of angles is between 100 to 180 degrees, this range for [C₃(mim)₂][NTf₂]₂ is between 70 to 140 degrees. In fact, in both of the studied GDILs, the obtuse angles are more probable. According to Fig. 8, the angles with the largest population for [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂ are about 120 and 140 degrees, respectively. It must be mentioned that, in [C₃(mim)₂][NTf₂]₂, although the angle at about 120 degrees is the angle with the largest population, the structures with the angle about 90 degrees are also important. Increasing this angle with increasing the linkage alkyl chain may be due to a balance between increasing the polarity alternation distance as a result of reduction of electrostatic interactions between the head groups and anions and flexibility of the linkage chain. The dependence of angle distribution on temperature is shown in Fig. S3 (ESI†) for [C₅(mim)₂][NTf₂]₂ in the temperature range of 298 to 333 K, which represents weak dependence on temperature. According to this figure, the maximum angle increases from 135 to 140, when the temperature increases from 298 to 333 K. This narrow increasing may be interpreted as a result of increasing the flexibility at higher temperatures.

Fig. 8 Distribution of the angles formed by two vectors pointing from the rings towards the center atom of the linkage chain for the studied GDILs at 333.15 K.

The spatial heterogeneity of the studied GDILs may be characterized by an uneven distribution of different sites in these systems and can be quantify by the heterogeneity order parameter (HOP) which is defined as^55,56

where N is the total number of sites in the system, r_ij is the distance between sites i and j corrected by the periodic boundary conditions, and σ = L/N^1/3 with L being the length of the cubic simulation box. HOP can be used to investigate the effect of linkage alkyl chain length between the rings on the spatial heterogeneity of the studied GDILs. Because a shorter distance between two sites has larger weight, the HOP value is larger when more sites are close in space. Fig. 9 shows the snapshots of the different sites of [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂ at 333.15 K and 1 atm. In this figure, the first row is for alkyl chain sites, the second row for anions, the third row for cation head sites and the last row for all atom sites. As this figure shows, all sites aggregate to form several spatially heterogeneous domains, but this aggregation is relatively more for the alkyl chains of the cations. The HOPs for alkyl chains, cation heads, and anions of the studied GDILs are given in Fig. 10. According to the study of Wang et al.,⁵⁶ for homogenously distributed ideal particles, the HOP approaches a constant value of 15.5 ± 0.2 after N ≥ 125 (see Fig. 4 of that paper). As Fig. 10 shows, although all of the alkyl chain, cation head and anion groups have HOP values lower than 15.5, and the HOPs for the linkage alkyl chain groups are much lower than those for the head groups and the anions, indicating that the alkyl chain groups distribute more heterogeneously compared to the cation heads and anions. Moreover, the HOP values decrease with the elongation of the linkage alkyl chain which shows that the [C₅(mim)₂][NTf₂]₂ system is more heterogeneous.

Fig. 9 Snapshots of [C₃(mim)₂][NTf₂]₂ (left column) and [C₅(mim)₂][NTf₂]₂ (right column) at 333.15 K and 1 atm. First row is for alkyl chain sites, second row for anions, third row for cation head sites and the last row for all atom sites.

Fig. 10 Heterogeneity order parameter (HOP) for alkyl chains, cation heads, and anions of the studied GDILs.

3.2.2. Transport properties. Besides the static equilibrium properties, the dynamical properties, such as self-diffusivity, ionic conductivity, and viscosity can be also calculated using MD simulations. Unlike the structural properties which have been discussed in the previous section, the dynamical properties of ILs depend on the time evolution of the system. In our previous work,²⁸ we measured and analyzed some transport properties of the [C_n(mim)₂][NTf₂]₂ GDILs including shear viscosity, diffusion coefficient, and electrical conductivity as a function of temperature and linkage alkyl chain length and then the other parameters such as transport numbers of cations and anions were calculated. In this section, to give a detailed description of a single particle motion in the studied GDILs from the molecular point of view, the mean-square displacement (MSD) of different groups and sites of these ILs has been calculated as a function of time. The plots of MSD were used to evaluate the self-diffusion coefficient of each ion of GDILs and the calculated diffusions can be used to estimate the transference numbers of cations and anions and also viscosities and electrical conductivities of the studied GDILs at 333.15 K and 1.0 atm.

For simulating the transport properties of these GDILs, the simulations were performed under 20% charge reduction on all atomic sites. To consider the polarization effects of highly polarizable ions, the partial charge of each atom must be reduced more than that calculated by methods such as CHELPG. In fact, using the overall charge as a fitting parameter in classical force fields describing ILs, as reduction of the ionic charges, will accelerate the dynamics and give a better agreement between classical force-field simulations and experimental results, especially for transport properties. In recent years, in simulations of IL properties, several authors used reduced charges by applying a scaling factor on all atomic sites of the IL.^33,43,57 It must be mentioned that using the nonpolarizable force field with intact partial charges (such as that carried out by Yeganegi et al.²⁶ on some DILs) underestimates the dynamics and overestimates the values of viscosities obtained by the simulation. In this case, just the trends of calculated diffusion coefficients and viscosities remain unchanged. More complete discussion in this subject has been given in our previous work on the simulation of linear tricationic ILs (LTILs).³⁴

The self-diffusion coefficient (D) of ion i can be evaluated from both the Green–Kubo formula⁵⁸ and also from the Einstein relation. The Einstein relation is the most commonly used method to calculate D. Indeed, if a diffusive regime can be established at relatively long simulation times, the values of diffusion coefficients can be calculated from the MSD plots. Each simulated MSD for the anions and cations was calculated from an 8 ns production run which was performed under 20% reduction of partial charges. The Einstein relation is as follows

where r_i(t) is the position of the center of mass of the molecule at time t. The self-diffusion coefficient can be estimated by plotting MSD as a function of time. To be sure that the simulation is long enough and the diffusive regime has been reached, the log(MSD) versus log(t) and also β versus t plots for the ions of the studied GDILs are shown in Fig. S4 (ESI†). The MSD of each ion generally has power law time dependence as

MSDαt^β (6)

The quantity β can be obtained aswhere β characterizes the type of motion presented in the system. For a correct estimation of the self-diffusion coefficient, the liquid has to be in a real diffusive regime. In the real diffusive regime, at very long times, β should reach a value equal to one, in which the MSD of the ion increases linearly with time. Accordingly, β values were obtained for the cations and anions of the studied GDILs at 333.15 and 1 atm from the log–log plots of the ions. The plots of β versus t are shown in Fig. S4 (ESI†) for the cations and anions of the studied GDILs to determine the diffusive regime. The results showed that, in each case, a truly diffusive regime has been reached after 3 ns except for the cation of [C₃(mim)₂][NTf₂]₂ in which the diffusive regime has been reached after 5.5 ns. The diffusion coefficients of ions of the studied GDILs were calculated using the MSD plots in the range of 3–8 ns (for [C₃(mim)₂][NTf₂]₂ in the range of 5.5–8 ns). Fig. 11 shows the MSDs of the cations and anions of the studied GDILs at 333.15 K and 1 atm. In addition, the values of the self-diffusion coefficients of the cations and anions of each GDIL are given in Table 3. A direct comparison between the simulated self-diffusion coefficients and the experimental values obtained by DOSY NMR spectroscopy in our previous work²⁸ is also given in this table. As this table shows, there is a relatively good agreement between the simulated and experimental values.

Fig. 11 Mean square displacements (MSDs) of the (a) cations and (b) anions of the studied GDILs at 333.15 K and 1 atm.

Table 3 The calculated values of the self-diffusion coefficients, transference numbers of cations and anions, and also viscosities and electrical conductivities of the studied GDILs at 333.15 K and 1.0 atm based on VACFs and also those calculated based on MSD plots. The values in parentheses are the relative percent deviations (Dev%)^a between the calculated and experimental values

Property	[C3(mim)2][NTf2]2			[C5(mim)2][NTf2]2
	Simu. (MSD)	Simu. (VACF)	Exp.	Simu. (MSD)	Simu. (VACF)	Exp.
a .
D + (10^−11 m^2 s^−1)	1.00 (8.93)	1.066 (2.91)	1.098	1.38 (−17.95)	1.202 (−2.73)	1.17
D − (10^−11 m^2 s^−1)	1.59 (11.67)	1.542 (14.33)	1.80	1.98 (0.50)	2.078 (−4.42)	1.99
% t+	55.71 (−1.53)	58.03 (−5.75)	54.87	58.23 (−7.97)	53.63 (0.56)	53.93
% t−	44.28 (1.86)	41.97 (7.00)	45.13	41.77 (9.33)	46.36 (−0.63)	46.07
Λ (mS cm^−1)	3.88 (10.00)	4.03 (6.51)	4.311	4.86 (−12.24)	4.46 (−3.00)	4.33
η (cP)	164.21 (−12.62)	165.40 (−13.44)	145.8	112.96 (−2.69)	114.00 (−3.64)	110.00

---

Furthermore, we calculated the self-diffusion coefficient (D) of the ions of the studied GDILs using velocity autocorrelation functions (VACFs). The VACF is defined using the following expression:

where the brackets represent the average over all times originating from t₀ within the trajectory, v_i represents the velocity of the center of mass of particle i, and t is the delay time of the correlation function. The VACFs for the cations and anions of the studied DILs at 333.15 and 1 atm are shown in Fig. S5 (ESI†). The calculated values of the self-diffusion coefficients based on VACF are listed in Table 3 and can be compared with the diffusions obtained from MSDs. However, the differences in the self-diffusion values obtained by these two methods are small and the values of the Einstein relation seem to be also acceptable.

Fig. 11 and also Table 3 show that for both studied GDILs, the anion diffuses faster than the corresponding dication. This is in line with the experimental results²⁸ due to the smaller size of the anion compared to the cation. The hydrodynamic radii of the anion/cation are 1.391/2.287 and 1.652/2.823 for [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂, respectively (see ref. 28). In addition, the self-diffusion of the cation and anions of [C₅(mim)₂][NTf₂]₂ is more than that of [C₃(mim)₂][NTf₂]₂ which is in agreement with the lower viscosity in the former. With an increase in the number of carbons in the linkage alkyl chain, the balance between two opposing effects, i.e. weakening the electrostatic attractions between the anion and cation as a result of the dispersion of the charge centers and also the enhancement of the vdW interactions, determines the dynamics of the ions. It seems that weakening the electrostatic attractions dominates and the diffusion of the [C₅(mim)₂][NTf₂]₂ increases.

A comparison among MSDs of the ring (the CR atom in the imidazolium cation), alkyl chain (the carbon atom in the middle of the alkyl chain), and anion (the nitrogen atom of the [NTf₂]⁻) for the studied GDILs is shown in Fig. S6 (ESI†) at 333.15 K and 1 atm. As this figure shows, for both studied GDILs, at large time scales, anions move faster than the atoms of the ring and these atoms, in turn, move faster than the atoms of the linkage alkyl chains. The time scales in which no crossing point has been observed are after 3 ns and 5 ns for [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂, respectively. This shows that the dynamics of these groups is identical at large time scales. Also, the MSDs show that the dynamics of the linkage alkyl chain and the ring in each of the studied GDILs are similar, which indicates that the motion of these groups is highly correlated with each other. However, there are some crossing points for these GDILs at short time scales. The crossing points are sensitive to the proper equilibration and trajectory averaging in the simulations for each IL.²⁶

Fig. S7 (ESI†) shows the MSDs of the carbon atoms of the linkage alkyl chain in the studied GDILs at 333.15 K and 1 atm. As shown in this figure, on moving from the rings toward the center of the linkage alkyl chain, the MSDs and thus self-diffusion coefficients decrease. In fact, the alkyl chain carbon atoms which are near the rings move faster than the carbons in the center of the linkage alkyl chain. It seems that the presence of the anions around the rings may be a possible reason for this result, due to the faster motion of the anions as discussed before. In other words, the anion tends to remain around the rings and the interaction between the anion and atoms of the chain decreases for the atoms that are far from the rings.

To investigate the temporal heterogeneity of the studied GDILs, the deviation of the self-part of the van Hove correlation function, G_s(r⃑,t),⁵⁹ from the Gaussian distribution of particle displacement and also the second-order non-Gaussian parameter, α₂(t),^60,61 were used. These quantities are defined as follow:

with

Although on macroscopic timescales the IL is homogeneous and there is no difference in the dynamical behavior of any two regions, in the simulation timescales (a few nanoseconds) there are significant differences in the dynamics of different regions.⁶²Fig. 12 shows the simulated non-Gaussian parameter associated with cation and anion displacement in the studied GDILs at 333.15 K and 1 atm. Larger values of α₂(t) can be related to the varying rates of both cationic and anionic motions, and their occurrence at longer times is due to the possibility of a non-diffusive mechanism (such as particle jump) even much after the onset of the diffusion timescale for a system.^63–65 As Fig. 12 shows, there is significant temporal heterogeneity for both the anion and cation in each of these GDILs, but this temporal heterogeneity is more for anions, which indicates that the anions have a higher degree of heterogeneity in their dynamics than cations. In addition, in the Fig. S4 (ESI†), the plots of β(t) = d[log(MSD)]/d[log(t)] versus t have been shown for the cations and anions of the studied GDILs to examine a diffusive regime. Fig. 12 represents that, although the maximum of α₂(t) occurs at shorter times (<1 ns) for the cation and anion of [C₅(mim)₂][NTf₂]₂ and the anion of [C₃(mim)₂][NTf₂]₂, the maximum of α₂(t) for the cation of [C₃(mim)₂][NTf₂]₂ occurs at a much longer time (>5 ns) which is in accordance with the results obtained in Fig. S4 (ESI†) in which β(t) become close to 1 at longer times to achieve the diffusive regime. However, as Pal and Biswas⁶³ mentioned, the presence of spatial heterogeneity (microscopic domain formation) is not a prerequisite for showing variation in motional rates, and temporal heterogeneity can exist even for much simpler systems.

Fig. 12 Non-Gaussian parameter associated with cation and anion displacement at 333.15 K and 1 atm for [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂.

The simulated self-part of the van Hove correlation function which is the probability distribution of displacement at time t, is shown in Fig. 13 for both of the ions in the studied GDILs at 333.15 K and 1 atm. Simulated functions correspond to the times at which the non-Gaussian parameters showed maxima. These times are mentioned in the accompanying legends in Fig. 13. In order to reflect the extent of deviation from Gaussian behavior of the particle displacements, the predicted Gaussian behavior is also shown in Fig. 13 at the same times. The Gaussian behavior can be accessed via the following relation:

Fig. 13 The simulated self-part of the van Hove correlation function for both of the ions in (a) [C₃(mim)₂][NTf₂]₂ and (b) [C₅(mim)₂][NTf₂]₂ at 333.15 K and 1 atm. Simulated functions correspond to the times at which the non-Gaussian parameters showed maxima. These times are mentioned in the accompanying legends. The dashed lines denote the calculations using Gaussian approximation.

As it can be seen from Fig. 13, G_s(r⃑,t) deviates from the Gaussian behavior for all of the studied ions. Also, the longer tail of the curve for the cation of [C₃(mim)₂][NTf₂]₂ shows relatively large amplitude displacements for the cation of this IL. This may be due to the shorter alkyl chain linkage between the rings in [C₃(mim)₂][NTf₂]₂.

Viscosity can be evaluated computationally from both the equilibrium MD (Green–Kubo) (GK) formalism, and Stokes–Einstein formula (it must be mentioned that there is also a non-equilibrium MD approach which is based on the response of the system to an external perturbation). In the Green Kubo method, the viscosity is calculated from the time integral over the stress autocorrelation function (SACF), i.e.. Generally, obtaining good statistics in the calculation of viscosity from the GK relation, requires averaging over long time trajectories to ensure that the integrand has decayed sufficiently to allow the upper limit on the integral to be truncated. In fact, the GK method is difficult to apply in simulations of highly viscous ILs since the correlation function converges slowly with time.^66,67 Thus, it seems that the estimation of the values of viscosities of the studied GDILs based on the Stokes–Einstein formula from the diffusion coefficient values can be also used as an approximation method. Viscosity can be estimated using the diffusion coefficients by the Stokes–Einstein equation

where D is the self-diffusion coefficient, r_h is the hydrodynamic radius, η is the viscosity, k_B is the Boltzmann constant, T is the absolute temperature and c is a constant. The value of c is different depending on the viscosity of the fluid. For large molecular size solutes in small molecular size solvent environments, c can be as high as 6. As the solute and solvent become more similar in size, especially in higher viscous liquids such as ILs, the value of c is reduced to about 4.³³ We used the value of c = 4 in all calculations. The effective radii, r_h, for the ions have been used from our previous work.²⁸ In this work, the viscosities of the studied GDILs have been calculated at 333.15 K and 1 atm using the Stokes–Einstein equation based on the diffusion coefficients calculated from both MSD and also VACF plots. The results are summarized in Table 3. As discussed in Section 3.1, the interaction energy between the cation and anions of [C₃(mim)₂][NTf₂]₂ is greater than that for [C₅(mim)₂][NTf₂]₂ which increases its viscosity. In fact, the interactions, including hydrogen bonds, electrostatic energies and vdW interactions determine the viscosity of an IL. Furthermore, we used the Green–Kubo method to calculate the viscosity values of both [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂ at 333.15 K and 1 atm (see Fig. S8, ESI†). The experimental viscosities of [C₃(mim)₂][NTf₂]₂ and [C₅(mim)₂][NTf₂]₂ at 333.15 K and 1 atm are 145.8 (cP) and 110.00 (cP). The calculated viscosities based on MSD, VACF plots (using Stokes–Einstein formula) and also the Green–Kubo method are 164.21, 165.40 and 166.19 cP (with the deviations −12.62%, −13.44% and −13.98%, respectively) for [C₃(mim)₂][NTf₂]₂. These values for [C₅(mim)₂][NTf₂]₂ are 112.96, 114.00 and 131.75 cP (with the deviations −2.69%, −3.64% and −19.77%, respectively). The results of both methods are acceptable, although the Stokes–Einstein formula just gives an approximation of the viscosities of these DILs.

The simulated self-diffusion coefficients can be used to calculate some other properties such as electrical conductivity, and also the transference numbers of cations and anions of the studied GDILs. The electrical conductivity was calculated from the Nernst–Einstein equation⁶⁸

where N denotes the number of ion pairs, e is the electronic charge, V is the volume of the system, k_B is the Boltzmann constant, T is the temperature and Z₋(Z₊) is the formal charge on the anion (cation). The electrical conductivity of the system depends on the diffusion of individual species. The applicability of the Nernst–Einstein equation for predicting the electrical conductivity of ILs has been investigated in some previous studies.^69–71 However, it must be mentioned that, due to the longer-range interactions and non-spherical ionic shapes of the studied GDILs, the results obtained by the Nernst–Einstein equation can be considered just as estimations.

The transference numbers were also estimated to study the contributions of the anions and cations to the transport of charges in these GDILs. The transference numbers show the relative mobility of the ions in an electrolyte. The transference numbers in a GDIL may be estimated from the diffusion coefficients of both cation and anion, according to

The results of electrical conductivities and also transference numbers of cations and anions in each of the studied GDILs are given and compared with corresponding experimental values in Table 3. As this table shows, there is relatively good agreement between calculated values in this work and those measured and obtained in our previous work.²⁸ Viscosity plays a main role in conductivity; the larger the viscosity, the lower the conductivity. The electrical conductivity of [C₅(mim)₂][NTf₂]₂ is greater than that of [C₃(mim)₂][NTf₂]₂ which is in contrast with the viscosity values of these ILs. Although, the greater ion size in [C₅(mim)₂][NTf₂]₂ can decrease the electrical conductivity, it seems that the viscosity is a dominant factor in the total electrical conductivity of these GDILs. The calculated values of cationic and anionic transference numbers of the studied GDILs are given in Table 3 at 333.15 K and at 1 atm. As this table shows, the transference number for the cation is greater than the anion which can be related to the larger diffusion coefficients of the cation. This result is in accordance with the results obtained in some previous studies for MILs,⁷² dicationic ILs²⁵ and tricationic ILs.³⁴ It must be mentioned that in electrolytes with t₊ > 0.5, the cations have a larger contribution to the ionic conductivity of the electrolyte than the anions.

4. Conclusions

The transport and interfacial properties and ionic characteristics of imidazolium-based short alkyl chain length DILs ([C_n(mim)₂][NTf₂]₂) as a function of linkage chain length and temperature have been studied experimentally in our previous studies. In this study, we obtained a fundamental understanding of the molecular basis of the macroscopic and microscopic properties of a bulk liquid phase of these GDILs using classical MD simulations and also quantum mechanical methods. Interaction energies, charge transfers and hydrogen bonds between the studied ion triplets were investigated by DFT calculations and also the AIM method. The results showed that the interaction energies between ion triplets in the studied GDILs are much greater than those of ion pairs in MILs which may be due to the stronger electrostatic attractions between the cation and the anions in GDILs and is responsible for their higher melting point and thermal stability. The results also showed that, in both of the studied GDILs, the interactions of oxygen atoms of the anions with the HA atom of each ring of the cation are stronger than the interactions of nitrogen atoms of the anions with the HA atoms. Also, it can be seen that the interactions of oxygen atoms of the anions with the HA atom of each ring are much greater than the interactions with the hydrogens of methylene groups connected to the rings.

The effect of linkage alkyl chain length on the thermodynamic, transport and structural properties of these GDILs have been investigated. The structural features of these GDILs were characterized by calculating the partial site–site RDFs and SDFs. In accordance with the results of DFT and AIM, RDFs show that, compared to the nitrogen atom, the oxygen atoms of the anions are more likely to form hydrogen bonds with the ring hydrogen atoms. The center of mass RDFs of cation–cation, anion–anion, and cation–anion for the studied GDILs show that to preserve the conditions of local electroneutrality, the ions of a given type prefer to be surrounded by ions of opposite charge. Alkyl chain–alkyl chain RDFs of the studied GDILs show the aggregation of the alkyl chains in [C₅(mim)₂][NTf₂]₂. In fact, with increasing the linkage alkyl chain length in [C₅(mim)₂][NTf₂]₂, the Coulombic interactions on the rings have less effect on the aggregation of the alkyl chain groups and the alkyl chain groups have more freedom to be closer to each other. In the case of [C₃(mim)₂][NTf₂]₂, the geometrical constraint restricts the aggregation of the alkyl chains due to the short length of these linkages. Also, the SDFs of these GDILs indicated that these ILs bend to maximize the electrostatic interactions of imidazolium rings with the anions and also to minimize the interactions between ions with the same charges.

The heterogeneity order parameter (HOP) has been used to describe the spatial structures of these GDILs and the distribution of the angles formed between two cation heads and the middle carbon atom of the linkage alkyl chain was analyzed in these ILs. According to the results, all sites aggregate to form several spatially heterogeneous domains, but this aggregation is relatively more for the alkyl chains of the cations. Also, for the studied GDILs, neither the straight chain model, nor the folded chain model exists and the nanoscale aggregates consist of a mixture of a majority of unfolded and a small portion of highly folded DILs.

To investigate the temporal heterogeneity of the studied GDILs, the deviation of the self-part of the van Hove correlation function, G_s(r⃑,t), from the Gaussian distribution of particle displacement and also the second-order non-Gaussian parameter, α₂(t), were used. G_s(r⃑,t) deviates from the Gaussian behavior for all of the studied ions. The longer tail of the curve for the cation of [C₃(mim)₂][NTf₂]₂ shows the relatively large amplitude displacements for the cation of this IL. This may be due to the shorter alkyl chain linkage between the rings in [C₃(mim)₂][NTf₂]₂. Although the maximum of α₂(t) occurs at shorter times (<1 ns) for the cation and anion of [C₅(mim)₂][NTf₂]₂ and the anion of [C₃(mim)₂][NTf₂]₂, the maximum of α₂(t) for the cation of [C₃(mim)₂][NTf₂]₂ occurs at a much longer time (>5 ns) which is in accordance with the results of β(t) plots in which β(t) for the cation of [C₃(mim)₂][NTf₂]₂ becomes close to 1 at longer times to achieve the diffusive regime.

MSD, self-diffusivity, viscosity, electrical conductivity and transference numbers of cations and anions have been computed for the studied GDILs at 333.15 K and at 1 atm. The simulated values were in good agreement with experimental data. The MSDs and self-diffusion values showed that weakening the electrostatic attractions dominates over the enhancement of the vdW interactions and therefore the diffusion of [C₅(mim)₂][NTf₂]₂ increases. The MSD plots showed that, for both studied GDILs, at large time scales, the anions move faster than the atoms of the ring and these atoms, in turn, move faster than the atoms of linkage alkyl chains. The results showed that, the cations have a larger contribution to the ionic conductivity of the electrolyte than the anions. Altogether, in this work, a better perspective of the structuring and dynamics of the studied GDILs has been obtained at a molecular level.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was supported financially by the Research Council of University of Isfahan. The authors are also indebted to Institute for Studies in Theoretical Physics and Mathematics (IPM) for providing us with the HPC Cluster to perform a part of the simulations.

References

1. A. Sanmartín Pensado, P. Malfreyt and A. l. A. Pádua, J. Phys. Chem. B, 2009, 113, 14708–14718.
2. B. Heggen, W. Zhao, F. Leroy, A. J. Dammers and F. Müller-Plathe, J. Phys. Chem. B, 2010, 114, 6954–6961.
3. T. G. Youngs and C. Hardacre, ChemPhysChem, 2008, 9, 1548–1558.
4. H. Sun, D. Zhang, C. Liu and C. Zhang, THEOCHEM, 2009, 900, 37–43.
5. B. Bhargava and M. L. Klein, J. Phys. Chem. B, 2011, 115, 10439–10446.
6. S. K. Chellappan, Y. Xiao, W. Tueckmantel, K. J. Kellar and A. P. Kozikowski, J. Med. Chem., 2006, 49, 2673–2676.
7. C.-M. Jin, C. Ye, B. S. Phillips, J. S. Zabinski, X. Liu, W. Liu and M. S. Jean'ne, J. Mater. Chem., 2006, 16, 1529–1535.
8. Z. Zeng, B. S. Phillips, J.-C. Xiao and J. N. M. Shreeve, Chem. Mater., 2008, 20, 2719–2726.
9. M. Palacio and B. Bhushan, J. Vac. Sci. Technol., A, 2009, 27, 986–995.
10. G. Yu, S. Yan, F. Zhou, X. Liu, W. Liu and Y. Liang, Tribol. Lett., 2007, 25, 197–205.
11. J. L. Anderson and D. W. Armstrong, Anal. Chem., 2005, 77, 6453–6462.
12. M. Qi and D. W. Armstrong, Anal. Bioanal. Chem., 2007, 388, 889–899.
13. K. Huang, X. Han, X. Zhang and D. W. Armstrong, Anal. Bioanal. Chem., 2007, 389, 2265–2275.
14. Z. Zhang, H. Zhou, L. Yang, K. Tachibana, K. Kamijima and J. Xu, Electrochim. Acta, 2008, 53, 4833–4838.
15. F. M. Nachtigall, Y. E. Corilo, C. C. Cassol, G. Ebeling, N. H. Morgon, J. Dupont and M. N. Eberlin, Angew. Chem., Int. Ed., 2008, 47, 151–154.
16. J. Y. Kim, T. H. Kim, D. Y. Kim, N.-G. Park and K.-D. Ahn, J. Power Sources, 2008, 175, 692–697.
17. C. Zafer, K. Ocakoglu, C. Ozsoy and S. Icli, Electrochim. Acta, 2009, 54, 5709–5714.
18. J. L. Anderson, R. Ding, A. Ellern and D. W. Armstrong, J. Am. Chem. Soc., 2005, 127, 593–604.
19. J. Pitawala, A. Matic, A. Martinelli, P. Jacobsson, V. Koch and F. Croce, J. Phys. Chem. B, 2009, 113, 10607–10610.
20. K. Ito, N. Nishina and H. Ohno, Electrochim. Acta, 2000, 45, 1295–1298.
21. B. Bhargava and M. L. Klein, J. Chem. Theory Comput., 2010, 6, 873–879.
22. D. Majhi and M. Sarkar, Phys. Chem. Chem. Phys., 2017, 19, 23194–23203.
23. N. S. Nooruddin, P. G. Wahlbeck and W. R. Carper, Tribol. Lett., 2009, 36, 147–156.
24. K. R. Lovelock, A. Deyko, J. A. Corfield, P. N. Gooden, P. Licence and R. G. Jones, ChemPhysChem, 2009, 10, 337–340.
25. S. Li, G. Feng, J. L. Bañuelos, G. Rother, P. F. Fulvio, S. Dai and P. T. Cummings, J. Phys. Chem. C, 2013, 117, 18251–18257.
26. S. Yeganegi, A. Soltanabadi and D. Farmanzadeh, J. Phys. Chem. B, 2012, 116, 11517–11526.
27. E. Bodo, M. Chiricotto and R. Caminiti, J. Phys. Chem. B, 2011, 115, 14341–14347.
28. M. Moosavi, F. Khashei, A. Sharifi and M. Mirzaei, Ind. Eng. Chem. Res., 2016, 55, 9087–9099.
29. M. Moosavi, F. Khashei, A. Sharifi and M. Mirzaei, J. Chem. Thermodyn., 2017, 107, 1–7.
30. M. Frisch, G. Trucks, H. Schlegel, G. Scuseria, M. Robb, J. Cheeseman and J. Montgomery, Gaussian 03, revision C. 02, 2008.
31. S. Tsuzuki, W. Shinoda, H. Saito, M. Mikami, H. Tokuda and M. Watanabe, J. Phys. Chem. B, 2009, 113, 10641–10649.
32. C. M. Breneman and K. B. Wiberg, J. Comput. Chem., 1990, 11, 361–373.
33. M. H. Ghatee, A. R. Zolghadr, F. Moosavi and Y. Ansari, J. Chem. Phys., 2012, 136, 124706.
34. E. Sedghamiz and M. Moosavi, J. Phys. Chem. B, 2017, 121, 1877–1892.
35. H. Kruse, L. Goerigk and S. Grimme, J. Org. Chem., 2012, 77, 10824–10834.
36. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
37. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354–360.
38. F. Biegler-Konig and J. Schonbohm, J. Comput. Chem., 2002, 23, 1489–1494.
39. W. Smith, T. Forester and I. Todorov, STFC, STFC Daresbury Laboratory, Daresbury, Warrington, Cheshire, WA4 4AD, UK, version, 2012.
40. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
41. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, 1989.
42. M. Salanne, Phys. Chem. Chem. Phys., 2015, 17, 14270–14279.
43. V. Chaban, Phys. Chem. Chem. Phys., 2011, 13, 16055–16062.
44. I. Leontyev and A. Stuchebrukhov, Phys. Chem. Chem. Phys., 2011, 13, 2613–2626.
45. W. Zhao, H. Eslami, W. L. Cavalcanti and F. Müller-Plathe, Z. Phys. Chem., 2007, 221, 1647–1662.
46. S. Aparicio and M. Atilhan, J. Phys. Chem. B, 2012, 116, 9171–9185.
47. K. Dong, S. Zhang, D. Wang and X. Yao, J. Phys. Chem. A, 2006, 110, 9775–9782.
48. P. A. Hunt, I. R. Gould and B. Kirchner, Aust. J. Chem., 2007, 60, 9–14.
49. M. Solimannejad, I. Alkorta and J. Elguero, Chem. Phys. Lett., 2007, 449, 23–27.
50. R. G. Parr and W. Yang, Annu. Rev. Phys. Chem., 1995, 46, 701–728.
51. M. Brehm and B. Kirchner, J. Chem. Inf. Model, 2011, 51, 2007–2023.
52. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
53. O. Russina, A. Triolo, L. Gontrani, R. Caminiti, D. Xiao, L. G. Hines Jr, R. A. Bartsch, E. L. Quitevis, N. Plechkova and K. R. Seddon, J. Phys.: Condens. Matter, 2009, 21, 424121.
54. A. Triolo, O. Russina, B. Fazio, R. Triolo and E. Di Cola, Chem. Phys. Lett., 2008, 457, 362–365.
55. Y. Ji, R. Shi, Y. Wang and G. Saielli, J. Phys. Chem. B, 2013, 117, 1104–1109.
56. Y. Wang and G. A. Voth, J. Phys. Chem. B, 2006, 110, 18601–18608.
57. B. Bhargava and S. Balasubramanian, J. Chem. Phys., 2007, 127, 114510.
58. J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids: With Applications to Soft Matter, Academic Press, 2013.
59. Z. Hu and C. J. Margulis, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 831–836.
60. J. Picalek and J. Kolafa, Mol. Simul., 2009, 35, 685–690.
61. A. Rahman, Phys. Rev., 1964, 136, A405.
62. S. S. Sarangi, W. Zhao, F. Müller-Plathe and S. Balasubramanian, ChemPhysChem, 2010, 11, 2001–2010.
63. T. Pal and R. Biswas, Theor. Chem. Acc., 2013, 132, 1348.
64. T. Pal and R. Biswas, J. Phys. Chem. B, 2015, 119, 15683–15695.
65. T. Pal and R. Biswas, J. Chem. Phys., 2014, 141, 104501.
66. N.-T. Van-Oanh, C. Houriez and B. Rousseau, Phys. Chem. Chem. Phys., 2010, 12, 930–936.
67. W. Jiang, T. Yan, Y. Wang and G. A. Voth, J. Phys. Chem. B, 2008, 112, 3121–3131.
68. J. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986.
69. S. U. Lee, J. Jung and Y.-K. Han, Chem. Phys. Lett., 2005, 406, 332–340.
70. C. Rey-Castro, A. Tormo and L. Vega, Fluid Phase Equilib., 2007, 256, 62–69.
71. C. Rey-Castro and L. F. Vega, J. Phys. Chem. B, 2006, 110, 14426–14435.
72. M. Kowsari, S. Alavi, M. Ashrafizaadeh and B. Najafi, J. Chem. Phys., 2008, 129, 224508.

---

Footnote

† Electronic supplementary information (ESI) available: The electrostatic charges for the atoms of the studied GDILs obtained at the B3LYP/6-31++G(d,p) level of theory by the CHELPG method, the simulated densities (in g cm⁻³) of the studied GDILs at different temperatures and at 1 atm and also some figures. See DOI: 10.1039/c7cp05681h

---