Tricationic Ionic Liquids: Structural and Dynamical Properties via

Article pubs.acs.org/JPCB

Tricationic Ionic Liquids: Structural and Dynamical Properties via Molecular Dynamics Simulations Elaheh Sedghamiz and Majid Moosavi* Department of Chemistry, University of Isfahan, Isfahan 81746-73441, Iran S Supporting Information *

ABSTRACT: Three imidazolium-based linear tricationic ionic liquids (LTILs) have been simulated to study their structural and dynamical properties and obtain a fundamental understanding of the molecular basis of the microscopic and macroscopic properties of their bulk liquid phase. The effects of temperature and alkyl chain length on the physiochemical, transport, and structural properties of these LTILs have been investigated. A nonpolarizable all-atom force field, which is a refined version of the Canongia Lopes and Paudua force field, was adopted for the simulations. Densities, mean square displacements, self-diffusivities, viscosities, electrical conductivities, and transference numbers have been presented for various ions from MD simulations. The detailed microscopic structures have been discussed in terms of radial distribution functions and spatial distribution functions. The results show that, similar to that in monocationic and dicationic ILs (MILs and DILs, respectively), the anions are mainly organized around the imidazolium rings. The diffusion coefficients of the studied LTILs are smaller than those of both MILs and DILs, with comparable viscosities. Unlike those of MILs and DILs, the diffusion coefficients of the cations and anions of the studied LTILs increase with an increase in the length of the alkyl chain between the rings for LTIL-1 and LTIL-2 but then decrease for LTIL-3, which is in a good agreement with the trend of viscosity data. The calculated transference numbers show that, similar to that in MILs and DILs, cations have a major role in carrying electric current in LTILs, but this role increases from MILs to LTILs.

1. INTRODUCTION In recent years, room-temperature ionic liquids (RTILs) have been placed among the most studied and attractive classes of new materials. ILs have wide range of applications, which are closely related to their unique properties, such as low vapor pressure,1 good electrical conductivity,2 electrical and chemical stability,3 wide liquid range,4 etc., making it and all-in-one material. Moreover, the most attractive feature of ILs is their tunability. Owing to the wide variety of cationic and anionic moieties, the physical and chemical properties of ILs can be considerably tuned for specific applications.5 Over the past few years, a large number of studies have focused on the correlation between the properties and structural changes of different groups of ILs.6−8 Most recently, it has been shown that simple modifications of traditional monocationic ILs (MILs) can produce multicationic ILs (especially dicationic ILs (DILs) and tricationic ILs (TILs)), which are more useful for a variety of applications. On the basis of the literature, it is evident that multifunctional ILs have a greater thermal stability and higher viscosities than those of the most common MILs, with almost the same solvent properties.9,10 Moreover, as the number of possible combinations of the cation and anion in multifunctional ILs is greater compared to that in MILs, it seems that their physicochemical properties, such as viscosity, density, thermal stability, melting point, and solubility, are more tunable.11−13 Recently, different © 2017 American Chemical Society

types of DILs (i.e., ammonium-based, phosphonium-based, and imidazolium-based ones) were synthesized and characterized by several groups.14−16 Anderson and co-workers11 were the first group that investigated the relationship between the structure and physicochemical properties of 39 DILs, including imidazolium- and pyridinium-based ILs. Because of their unique characteristics, DILs have rapidly found some special applications, for example, as a stationary phase for gas chromatography,17 as a solvent for high-temperature organic reactions,18 as high-temperature lubricants,19 and as electrolytes in secondary batteries20 and dye-sensitized solar cells.21 As a special class of multicationic ILs, TILs comprise three charge centers with three singly charged cations linked to each other by two spacers, which are usually alkyl chains. In their crystalline phase, a singly charged anion is associated with each positive charge center. TILs have a trigonal or linear geometry. They are thermally highly stable and their liquid temperature range is more than about 300 °C.9,13 Linear TILs (LTILs) have some advantages over the trigonal type. For example, all of them are RTILs due to their higher structural flexibilities compared to those of the trigonal ones. This structural flexibility would give the best tunability to the LTILs in Received: October 25, 2016 Revised: February 1, 2017 Published: February 2, 2017 1877

DOI: 10.1021/acs.jpcb.6b10766 J. Phys. Chem. B 2017, 121, 1877−1892

Article

The Journal of Physical Chemistry B

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

terms of their physiochemical properties.13 In spite of the limited number of synthesized TILs, some of their applications have been explored. They have been used in all areas of separation science, including extractions, gas chromatography, and supported liquid membranes.22−24 In addition, they have shown great performance as ion-pairing reagents for the ultratrace detection of anions25 and in electro wetting applications.9 Molecular dynamics (MD) simulation has been an important tool in the field of IL research for the interpretation of

experimental observations and also for designing new IL systems for specific applications. To date, most theoretical research in the field of ILs is related to MILs and to some extent to DILs. Bodo and Caminiti26 reported the molecular mechanics and ab initio calculations of geminal di-imidazolium bis (trifluoromethylsulfonyl)amide ILs in the gas phase. They investigated the preferred geometries of the ionic complexes and reported the bending of the alkyl chain in DILs for maximizing the electrostatic interactions between cations and anions. In their second paper, Bodo et al.27 performed MD 1878


Article


propylimidazolium)imidazolium (LTIL-1), 1-(1′-methyl-3′hexylimidazolium)-3-(1″-methyl-3″-hexylimidazolium)imidazolium (LTIL-2), and 1-(1′-methyl-3′-decylimidazolium)3-(1″-methyl-3″-decylimidazolium)imidazolium (LTIL-3) were obtained as ±2.861e, ±2.5646e, and ±2.686e, respectively (which means the net charge on each anion would be approximately −0.95e, −0.85e, and −0.89e, respectively). The partial charges on each atom of the studied LTILs, calculated by the CHELPG method, have been given in Table S1 of the Supporting Information. 2.2. Force-Field Parameters. In recent years, several force fields have been developed for atomistic description of different groups of ILs. Among them, the nonpolarizable all-atom force field developed by Canongia Lopes and Paudua (CL&P)33−36 is the most widely used force field. However, the ability of nonpolarizable force fields to accurately predict all physical properties of ILs is less established. They greatly underestimate the ionic dynamics, but at the same time, they can reproduce more accurate thermodynamic and structural properties within a 3% deviation.37 In contrast, many-body polarizable force fields are more successful in predicting the dynamic properties of ILs, but they have a noticeably higher computational cost compared with nonpolarizable force fields.38 In the case of all-atom force fields in which the computed dynamic properties of ILs converge slowly, the use of polarizable force fields for simulation of ILs is not feasible except for ILs with limited number of atoms. In this situation, to improve the nonpolarizable force fields, we can 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.37,39 Schmidt et al.40 used the Blochl method, which is a reciprocal space approach for the evaluation of charges, to assess the behavior of the ions in the liquid state directly, unlike other methods. They concluded that just on the one-ion-pair level the Blochl method yields results that agree well with those of the standard procedures, such as CHELPG. In fact, their analysis supported the idea of 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. Therefore, 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 recent years, in simulations of IL properties, several authors used reduced charges by applying a scaling factor on all atomic sites of the IL. For example, Ghatee et al.41 reduced the calculated partial charges by 20% for [Cnmim][I], with n = 4, 6, 8 using the CHELPG method, leading to a faster dynamics. (Their calculated net charges using the CHELPG method were ±0.8750e, ±0.8707e, and ±0.8671e, respectively). Bhargava and Klein42 have also used the same procedure. In addition, Chaban37 studied the effect of polarizability on different properties of two MILs using different scaling factors for reducing the partial charges on all atomic sites. It must be mentioned that using the nonpolarizable force field with intact partial charges (such as that done by Yeganegi et al.31 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

simulations of the bulk liquid phase of geminal di-imidazolium bis(trifluoromethylsulfonyl) imide ILs. They reported longrange structural ordering of DILs and compared the bulk liquid structure with both gas-phase structures and the structures of the monoimidazolic counterparts. Vitorino et al.28 studied the vaporization of DIL [C3(mim)2] [NTf2]2 by means of the ioncyclotron resonance mass spectroscopy technique and MD simulation. Their results showed that the enthalpy of vaporization of DILs is not as high as expected and is comparable to that of MILs. Bhargava and Klein29 studied the formation and organization of micellar aggregates in an aqueous solution of DIL 1,3-bis(3-decylimidazolium-1-yl) propane bromide. In their study, the structure of the vapor−liquid interface of the solution and micellar aggregates and also the distribution of the counterions in water were discussed. In another paper, Bhargava and Klein30 studied the structure and organization of the ions in aqueous solutions of the same IL by coarse-grained MD simulations. Their results show that unlike the aggregates formed by MILs, which are spherical, those of DILs display a near-hexagonal arrangement. Yeganegi et al.31 investigated the structure and dynamics of nine geminal DILs. They studied the effects of the alkyl chain length of the cation and also the type of anion on the diffusion coefficients. Their results indicated that, unlike those of MILs, the anions and cations of DILs distribute homogeneously and that anions have a major role in the transfer of electric current. As the authors are aware, there is no available computational research on the properties of TILs. There is no doubt that deeper and broader knowledge of the physiochemical properties of TILs at a molecular level is essential for finding their unique features and extending their applications to different areas and also improving their rational design. MD simulations can provide valuable details on the structure, dynamics, and solvation properties of these new types of ILs. In this work, for the first time, three imidazolium-based LTILs with different alkyl chains between the rings (see Figure 1) have been simulated to obtain a fundamental understanding of the molecular basis of the microscopic and macroscopic properties of their bulk liquid phases. In addition, the effect of alkyl chain length on the physiochemical properties of these LTILs has been investigated. In Section 2, details of the computational methods, including ab initio calculations, forcefield parameters, and details of the MD simulations, are presented. In Section 3, results of bulk thermodynamic, transport, and structural properties are presented and have been discussed in detail.

2. COMPUTATIONAL METHODS 2.1. Ab Initio Calculations. The structures of the isolated cations 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.32 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 cation and anion charges were set at +3 and −1, respectively. Thereafter, the structure of each LTIL (cations and anions together) was optimized using B3LYP/631*G. Schematic structures of the studied LTILs are shown in Figure 1. 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 net charges on the cations and anions of 1-(1′methyl-3′-propylmidazolium)-3-(1″-methyl-3″1879


Article


Table 1. Comparison between the Simulated Densities of the Studied LTILs Using Original Charges and Charges Reduced by 20% with Both Experimental Values13 and Those Predicted Based on GCM49 at 296 K and 1 atm % dev. (with original charges)

ρ (g cm−3) LTIL-1 LTIL-2 LTIL-3

exp.

GC

simu. (with original charges)

simu. (with 20% charge reduction)

from exp.

from GC

from exp.

from GC

1.54 1.57 1.65

1.62 1.490 1.405

1.64 1.48 1.43

1.618 1.45 1.41

−6.49 5.73 13.33

−1.23 0.67 −1.78

−5.06 7.64 14.55

0.12 2.68 −0.36

ature of the systems was increased at intervals of 100 K up to 650 K. At each temperature, the simulation was performed for 500 ps under constant NPT conditions. Then, the temperature was decreased to 450 K at steps of 100 K. The simulations were continued at different temperatures, from 450 to 296 K. At each temperature, 500 ps was allocated for system equilibration under constant NVT with the Nosé−Hoover thermostat, followed by 500 ps under constant NPT conditions. The starting point of each NPT simulation was an equilibrated configuration of the NVT simulation. After the heating/cooling procedure, a 2 ns NPT simulation was carried out for data collection. Finally, the lengths of the sides of the simulation boxes were 52.64, 55.80, and 58.09 Å for LTIL-1, LTIL-2, and LTIL-3, respectively. The mean square displacements (MSDs) at 303 K were obtained by a 1 ns simulation under NPT conditions, followed by longer production times to reach the diffusive region.

remain unchanged. Thus, for simulating the transport properties of LTILs, the simulations were performed under 20% charge reduction on all atomic sites, which means that the total net charge on each anion of LTIL-1, LTIL-2, and LTIL-3 would be −0.76e, −0.68e, and −0.71e, respectively. In this work, a refined force field for imidazolium-based ILs, which was basically developed by Tsuzuki et al.,43 was employed to simulate the LTILs, namely, LTIL-1, LTIL-2, and LTIL-3. This force field is a refined version of the CL&P force field for imidazolium-based ILs, obtained on the basis of ab initio calculations. The potential energy of the systems was modeled using the following standard functional form Vtot =

∑

kb(r − req)2 +

bonds

∑

kθ(θ − θeq)2

angles 3

+

∑ ∑ dihedrals i = 1

Vi (1 + cos(iφ + fi )) 2

⎧ ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ qiqj ⎫ σij ⎥ ⎪ ⎪ ⎢ σij ⎬ ∑ ⎨4εij⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥ + r ⎝ rij ⎠ ⎦ 4πε0rij ⎪ j>1 ⎪ ⎭ ⎩ ⎣⎝ ij ⎠

3. RESULTS AND DISCUSSION 3.1. Temperature Dependence of Density. Density is one of the most fundamental properties of ILs and also one of the most accurate forms of experimental data. Because of its simplicity of calculation and the availability of the very accurate experimental values, density is usually the most common property used for validating the force-field parameters in MILs and DILs.47 Unlike that for MILs, the experimental density data for TILs is scarce, and only one set of data can be used to validate the proposed force field. In this regard, the group contribution method (GCM), which is developed to predict the thermodynamic properties of a large set of ILs, seems to be an alternative source in cases in which sufficient or accurate dataset are not available.48,49 The densities of three studied LTILs calculated by NPT simulations at 296 K and 1 atm have been reported and compared to the experimental data of Wanigasekara et al.13 in Table 1 and Figure 2. The deviations are rather high, although none of the force-field parameters have been adjusted to match the experimental data. As shown in Table 1 and Figure 2, the experimental densities of the LTILs increase with an increase in the number of carbon atoms. The results of our simulations show an opposite trend to these reported experimental data but are compatible with the trend resulting from GCM.49 Previous experimental and simulation studies on both MILs50−52 and DILs10,31 show that the higher the alkyl chain length the lower the density. It seems that the values reported by Wanigasekara et al.13 are not reliable, and this may be attributed to the impurities, including the water content, in the samples or the method of measurement in their work. As shown in Table 1, the densities resulting from simulation are in good agreement with the densities predicted by GCM49 within statistical error of the simulated values. The temperature dependence of density for the studied LTILs have been given in Figure 3a. As shown in this figure, the simulated densities of LTILs decrease

N−1 N

+

∑ i=1

% dev. (with 20% charge reduction)

(1)

The potential energy associated with the bonds and angles is described in harmonic terms, the dihedral torsion energy is represented by a series of cosines, and nonbonded interactions are given by the Lennard-Jones sites and Coulombic interactions (calculated using the Ewald summation method) between partial point charges placed on the atomic sites. The electrostatic charges, bond lengths, bond angles, and dihedrals for the cationic part were taken from our ab initio calculations. For the [NTf2]− anion, we have used electrostatic charges from our ab initio calculations, and all other parameters were taken from the force field.43 2.3. Simulation Details. An ion quartet (one cation and three anions) of each LTIL, with a geometrical structure optimized via ab initio calculations, was replicated to obtain a simulation box containing 125 ion quartets (125 cations and 375 anions), with a total number of 11 875, 14 125, and 17 125 atoms for LTIL-1, LTIL-2, and LTIL-3, respectively. For each IL system, the simulation was started at 350 K at ambient pressure (1.01325 × 105 Pa), using the DL_POLY program, version 2.17.44 Each simulation was done with a time step of 1 × 10−4 ps for the initial 400 ps and then followed by a time step of 1 × 10−3 ps. The Nosé−Hoover thermostat and barostat,45 with 0.1 and 2 ps relaxation times, respectively, were used to control the temperature and pressure. Long-range Coulombic interactions were handled using the Ewald summation46 method, with a precision of 1 × 10−5. A cutoff distance of 25.0 Å was chosen for the Lennard-Jones interactions. The equations of motion were solved using the Verlet−Leapfrog integration algorithm,46 and cubic periodic boundary conditions were applied. The simulations were continued until the systems reached an initial equilibrium. Thereafter, the temper1880


Article


linearly with increasing temperature in the range of 296 to 356 K. The decrease in the density of LTILs with increasing temperature is similar to that for traditional MILs50−53 and also DILs,10,31 but the densities of LTILs are higher than those of MILs with comparable alkyl chain lengths. It is usual to use nonpolarizable force fields in simulating different properties of ILs and improve the results by reducing the partial charges on all atomic sites. For example, Ghatee et al.41 showed that reducing the atomic charges by 20% of the original charges improves the simulated densities in a series of imidazolium-based iodide MILs. Chaban and Prezhdo54 showed that using a scaling factor of 0.8 (i.e., reducing the charges by 20%) gives the simulated densities very close to experimental values for two different MILs, that is, 1,3dimethylimidazolium bis(trifluoromethane sulfonyl)imide ([MMIM][TFSI]) and N-butylpyridinium bis(trifluoromethane sulfonyl)imide ([BPY][TFSI]). Yeganegi et al.31 used the all-atom force field of Lopes et al.35 without any modification for a series of DILs, that is, Cn(mim)2X2, with n = 3, 6, 9 and X = PF6−, BF4−, Br−. They just made a small modification in the cation atomic charges for the short alkyl chain compounds. In this work, we try to examine the effect of reducing charges by using a scaling factor. As shown in Table 1, although reducing the charges by 20% results in better agreement between the simulated densities and the GCM values for LTIL-1 and LTIL-3, it deteriorates the agreement with LTIL-2. Altogether, the simulated densities with the original charges are acceptable for all of the studied LTILs, and the deviation is less than 2%. Long-range interactions can lead to significant dependencies on the system size in computer simulations of condensed matter under periodic boundary conditions.55 This finite-sized effect could be even more significant in simulations of ILs. There exist a large number of studies on the different properties of ILs, which used different system sizes (from 96 to 800 ILs), depending on the type of IL, the used software, and also the applied force field.33,37,56−60 From these studies, it can be concluded that 125 ILs in the simulation box could give reasonable results, and the simulation results are not subject to finite-size effects. However, because of the novelty of simulation of LTILs, we tried to investigate a larger system with 216 ILs (20 520 atoms in a cubic box) for LTIL-1. Table 2 shows the effect of system size (the number of ion quartets) on the simulated densities of LTIL-1 using original charges and charges reduced by 20% with both experimental values13 and those predicted based on GCM49 at 296 K and 1 atm. The results show that increasing the system size is not significant and the densities obtained from 125 ILs seem to be reliable. Thermal expansion coefficients can be evaluated using the variation of density with temperature. These coefficients for the studied LTILs at a constant pressure were easily obtained from linear fits of the density data using eq 2

Figure 2. Comparison between the simulated densities and experimental values13 and also the densities predicted by GCM49 for the studied LTILs at 296 K and 1 atm. Lines are drawn to guide the eyes.

ρ (g cm−3) = a + bT αP (K−1) = −(1/ρ)(∂ρ /∂T )P = −b/(a + bT )

(2)

−1

where αP is the thermal expansion coefficient in K and a and b are the fitting parameters of eq 2. The variation of αP for the studied LTILs has been illustrated in Figure 3b. To compare the thermal expansion coefficients of LTILs with those of MILs and DILs, the values of αP for MIL [C3mim][NTf2]61 and DIL [C3(mim)2][NTf2]262 have been also shown in Figure 3b. According to this figure, the values of αP for the studied LTILs

Figure 3. (a) Temperature dependence of the simulated densities and (b) calculated isobaric expansion coefficients, αP, for the studied LTILs. The isobaric expansion coefficients of one MIL ([C3mim][NTf2])61 and one DIL (C3[mimNTf2]2)62 have been given for comparison. 1881


Article


Table 2. Effect of System Size (the Number of Ion Quartets) on the Simulated Densities of LTIL-1 Using Original Charges and Charges Reduced by 20% with Both Experimental Values13 and Those Predicted Based on GCM49 at 296 K and 1 atm % dev. (with original charges)

ρ (g cm−3)

% dev. (with 20% charge reduction)

system size (number of ion quartets)

exp.

GC

simu. (with original charges)

simu. (with 20% charge reduction)

from exp.

from GC

from exp.

from GC

125 216

1.54 1.54

1.62 1.62

1.640 1.636

1.618 1.609

−6.49 −6.23

−1.23 −0.98

−5.06 −4.48

0.12 0.68

Table 3. Calculated Diffusion Coefficients, Transference Numbers, and Hydrodynamic Radii of Cations and Anions and Also Viscosities and Electrical Conductivities of the Studied LTILs at 303 K and 1.0 atm LTILs

D+ (10−13 m2 s−1)

D− (10−13 m2 s−1)

t+

t−

Λ (10−5 S m2 mol−1)

r+ (Å)

r− (Å)

ηexp. (cSt)

ηsimu. (cSt)

LTIL-1 LTIL-2 LTIL-3

5.74 19.15 4.88

4.00 17.50 3.70

0.81 0.77 0.80

0.19 0.23 0.20

2.35 8.3 2.03

5.69 6.01 6.57

3.33 3.33 3.33

1200 372 1800

1004.7 284.9 1255.3

are lower than those for the considered MIL and DIL, which shows that extension of the volume of an LTIL is lower than that for both the MIL and DIL, in fact, the order of αP is as follows: MIL > DIL > TIL. Also, all of the ILs have less expansion capacity compared to that of traditional organic solvents.63 As can be seen in Figure 3b, the thermal expansion coefficients of these LTILs increase with increasing alkyl chain length of the cation and also temperature. Increasing the alkyl chain length causes dispersion of the positive charge center of the cation and weakening of the electrostatic attraction between the cations and anions, and this may be a reason for the increasing αP with increasing alkyl chain length in LTILs. 3.2. Transport Properties: Diffusion, Viscosity, and Electrical Conductivity. According to the previous studies on different groups of ILs, it has been shown that the experimental measurements of the transport properties of these compounds suffer some issues, such as impurity and water content of the samples.64 On the other hand, several factors are important in the quality of the results of the simulation of the transport properties of ILs, including the need to perform large-scale simulations, desired large simulation times, using suitable forcefield parameters, and, perhaps the most important of all, correction of the charges in the force fields when we use nonpolarizable ones.7 Nonpolarizable force fields overestimate the values of viscosities obtained by simulation. It is shown that reducing the atomic charges affects the long-range interactions and thus increases the diffusion coefficient and ion mobility.54 The MSD can be calculated via the following relation MSD =

N ⎞ 1⎛ ⎜⎜∑ |ri(t ) − ri(0)|2 ⎟⎟ = ⟨|ri(t ) − ri(0)|2 ⟩ N ⎝ i=1 ⎠

β (t ) =

1 d lim ⟨|ri(t ) − ri(0)|2 ⟩ 6 t →∞ dt

(5)

Obviously, when β approaches unity, the system exhibits diffusive motion. Figure S1 shows the evolution of β with time for the middle ring, side rings, alkyl chain, and anion for LTIL1, LTIL-2, and LTIL-3 at 303 K and 1 atm. As shown in this figure, the systems appear to achieve a diffusive region in the range of 10−30 ns. The evolution of the β parameter with time for the cation and anion of LTIL-1 at different temperatures up to 8 ns has been shown in Figure S2. According to this figure, diffusive motion is achieved much more rapidly at higher temperatures. Thus, for reliable estimation of the self-diffusion coefficients at low temperatures, longer simulation times are needed. The diffusion coefficients of the anions and cations of the studied LTILs, which were calculated using the slopes of the lines fitted to the MSDs in the range of 10−30 ns, are presented in Table 3 at 303 K and 1 atm. The MSDs of the CA atom in the middle ring of the imidazolium cation and the nitrogen atom of the [NTf2]− anion for each of the LTILs have been given in Figure 4. According to this figure, all MSDs are linear in the range of 10−30 ns. The calculated diffusion coefficients of the cations and anions in the ranges of 10−15, 14−20, 19−25, 24−30, and 10−30 ns for the studied LTILs have been given in Table 4. As shown in this table, there is good agreement between the simulated diffusion coefficients over different time ranges. Yeganegi et al.31 showed that the diffusion coefficient of a DIL is smaller than that of a MIL. From the values of diffusion coefficients in Table 3, it is clear that the diffusion coefficient of a LTIL is smaller than that of both MILs and DILs, with comparable viscosities. A number of factors can influence the diffusion coefficients of the anions and cations in ILs. Some of these factors are the relative cationic and anionic sizes; geometrical shape; mass; strength of local van der Waals interactions between the cations and anions; strength of Coulombic interactions between the cation and anion, which depend on the localization of charge on the ions; and the surface electrical charge density on the ions.65 Yeganegi et al.31 showed that contrary to the observed trend in MILs, the diffusion coefficients of DILs increase with increasing alkyl chain length due to a balance between two opposing effects of chain lengthening on the dynamics of ions. It seems that a balance between two opposite factors determines the properties of ILs. Tokuda et al.6 showed that the addition of a −CH2− unit to the cation of an IL decreases

(3)

where ri is the vector coordinate of the center of mass of ion, i. MSDs of the anion and cation were obtained by performing simulations under 20% charge reduction on all atomic sites. Self-diffusion coefficients of the ions, Di, have been obtained from the Einstein relation through the slopes of the MSDs in the long time regime via the eq 4 Di =

d log⟨Δr 2⟩ d log(t )

(4)

It is of critical importance that the system exhibits true diffusive motion. To test for diffusive motion, parameter β is computed over a range of time scales, where β is defined according to 1882


Article


et al.,13 who also extensively discussed the effect of alkyl chain length on the viscosity of these LTILs. To explore the dynamics of the middle and side rings and the alkyl chain atoms of the cations of the studied LTILs, the MSDs of the middle ring (the CA atom in the middle ring of the imidazolium cation), the side ring (the CY atom in the side ring of the imidazolium cation), the alkyl chain (the carbon atom in the middle of the alkyl chain (C2 for LTIL-1, C4 for LTIL-2, and C6 for LTIL-3)), and the anion (the nitrogen atom of [NTf2]−) for each of the studied LTILs at 303 K and 1 atm have been compared in Figure 5. Table 5 shows the values of the calculated diffusion coefficients for the anion, middle ring (the CA atom in Figure 1), side ring (the CY atom in Figure 1), and alkyl chain (the carbon atom in the middle of the alkyl chain) from the slopes of the calculated MSDs in Figure 5 for the studied LTILs at 303 K and 1.0 atm. As shown in Figure 5, there is no crossing point in the MSDs of the middle and side rings and the alkyl chain atoms of the cations and also the MSDs of the anions for each of the studied LTILs, and these MSDs are almost parallel to each other. This shows that the dynamics of these groups is identical. According to Figure 5 and also Table 5, the MSDs and diffusion coefficients of the alkyl chain, middle ring, and side ring in each of the studied LTILs are similar, which indicates that the motion of these groups is highly correlated to each other. Viscosity can be evaluated computationally from both the equilibrium MD (Green−Kubo (GK)) formalism and the Stokes−Einstein formula (it must be mentioned that there is also a nonequilibrium MD approach, which is based on the response of the system to an external perturbation). In the GK method, the viscosity is calculated from the time integral over the stress autocorrelation function, that is, ∞ 1 η = Vk T ∫ ⟨Pαβ(t )Pαβ(0)⟩dt . Generally, obtaining good staB

Figure 4. MSDs of (a) the CA atom in the middle ring of the imidazolium cation and (b) the nitrogen atom of the [NTf2]− anion for all studied LTILs at 303 K and 1 atm.

Table 4. Calculated Diffusion Coefficients of the Cations and Anions from the Slopes of the MSD Plots in Figure 4 in the 10−15, 14−20, 19−25, 24−30, and 10−30 ns Time Ranges for the Studied LTILs at 303 K and 1.0 atm LTILs molecules

10−15 ns 14−20 ns 19−25 ns 24−30 ns 10−30 ns

D− (10−13 m2 s−1) LTIL-1 3.86 LTIL-2 15.54 LTIL-3 3.83 D+ (10−13 m2 s−1) LTIL-1 5.45 LTIL-2 22.20 LTIL-3 4.76

4.03 19.02 3.79

4.07 16.64 3.69

4.02

5.80 18.07 5.03

5.66 20.80 5.33

5.72

3.62

4.66

0

tistics in the calculation of viscosity from the GK relation requires averaging over very 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, as the correlation function converges very slowly with time66,67 VanOanh et al.66 studied the effect of simulation time on the convergence of the simulation of [emim][NTf2] in the GK method. Their results showed that the viscosity converges only after about 30 ns at a high temperature (at 500 K, with a viscosity of 3.57 mPa s). However, they showed that it is much more difficult to reach the convergence of viscosity as the temperature decreases because the dynamics slows with decreasing temperature. They showed that at 293 K a simulation time of more than 100 ns is required. In the case of the studied LTILs, which have much higher viscosities than those of [emim][NTf2], it is obvious that longer simulation times are required (longer than 100 ns). Thus, in this work, the Stokes−Einstein formula, which only leads to an estimation of the viscosity, has been used. The Stokes−Einstein relation is as follow

4.00 17.50 3.70 5.74 19.15 4.88

the electrostatic attraction between the cation and anion. On the other hand, an increase in −CH2− units enhances the van der Waals interactions by means of alkyl chain−ion inductive forces and hydrocarbon−hydrocarbon interactions. In the case of the studied LTILs, it can be observed that the diffusion coefficients of the cations and anions increase with an increase in the length of the alkyl chain from LTIL-1 to LTIL-2 but then decrease for LTIL-3. This is in excellent agreement with the experimental trends of viscosity data observed by Wanigasekara

η=

kBT cπDrh

(6)

where η is the viscosity, D is the self-diffusion coefficient, rh is the hydrodynamic radius, kB 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 solutes of large molecular size in solvent environments of small molecular 1883


Article


Table 5. Calculated Diffusion Coefficients for the Anion, Middle Ring (the CA Atom in Figure 1), Side Ring (the CY Atom in Figure 1), and Alkyl Chain (the Carbon Atom in the Middle of the Alkyl Chain (C2 for LTIL-1, C4 for LTIL-2, and C6 for LTIL-3) from the Slopes of the Calculated MSDs in Figure 5 for the Studied LTILs at 303 K and 1.0 atm D (10−13 m2 s−1) LTIL molecules LTIL-1 LTIL-2 LTIL-3

anion

D

5.70 19.95 4.00

alkyl chain

D

4.95 21.03 3.78

Dmiddle ring

Dside rings

3.90 17.30 3.65

5.56 23.86 3.68

sphere that diffuses at the same rate as a solute.68 To calculate the viscosities of the studied LTILs, the hydrodynamic radii of these ILs must be estimated. To estimate the hydrodynamic radius of these compounds, the method proposed by MacFarlane et al.2 was used. According to this method, the molecular volumes were computed using DFT at the B3LYP/6311++G(d,p) level of theory. Normally, within this approach the volumes are computed with an accuracy of 10%. The effective ionic radii are calculated, assuming that the ions studied have a spherical shape 4 3 V theory = πR ion (7) 3 It should be noted that for most of the ions, the effective radius defined by this method is an underestimation of their true hydrodynamic radius. The hydrodynamic radii of the cations and anions of the studied LTILs have been given in Table 3. As shown in this table, the hydrodynamic radii of the cations increase with increasing length of the alkyl chain in these LTILs. As expected, the hydrodynamic radius of the cation is greater than that of the anion. The calculated viscosities of the studied LTILs at 303 K and 1 atm show good agreement with the experimental values.13 Molar electrical conductivity (Λ) can be evaluated from the Nernst−Einstein equation, using the diffusion coefficients of the cation and anion of the IL, as follow Λ=

NAe 2(V +Z+2D+ + V −Z −2D−) kT

(8) −1

where Λ is the molar electrical conductivity (S cm mol ), NA is Avogadro’s number, e is the electron charge (1.602 × 10−19 coulombs), V+ and V− are the number of ions per formula unit, Z+ and Z− are the charges of the ions, D+ and D− are the diffusion coefficients of the cation and anion (cm2 s−1), respectively, k is the Boltzmann constant, and T is the absolute temperature. As expected, the trend of the molar electrical conductivities of the studied LTILs is in contrast to that of the viscosity values of these ILs. LTIL-2, which has the lowest viscosity, has a greater molar electrical conductivity. Ion transport number or transference number is the fraction of the total current carried in an electrolyte by a given ion. Differences in transport numbers arise from differences in electrical mobility. The transport numbers of cations and anions of the studied LTILs were evaluated using their diffusion coefficients as 2

Figure 5. Comparison of the MSDs of the middle ring (the CA atom in the middle ring of the imidazolium cation), side ring (the CY atom in the side rings of the imidazolium cation), alkyl chain (the carbon atom in the middle of the alkyl chain (C2 for LTIL-1, C4 for LTIL-2, and C6 for LTIL-3)), and anion (the nitrogen atom of [NTf2]−) for (a) LTIL-1, (b) LTIL-2, and (c) LTIL-3 at 303 K and 1 atm.

size, c can be as high as 6. As the solute and solvent become more similar in size, especially in liquids of higher viscosity, such as ILs, the value of c is reduced to about 4.65 We used the value of c = 4 in all calculations. The hydrodynamic radius can refer to the Stokes−Einstein radius of a solute, which is defined as the radius of a hard

t+ = 1884

3D+ D− − , t− = + 3D + D 3D + D− +

(9)


Article

The Journal of Physical Chemistry B Transport numbers of the cations and anions in each of the studied LTILs have been given in Table 3. As shown in this table, the cation transport number is larger than the anion transport number, which is expected, according to the trends of diffusion coefficients in these ILs. The calculated transport numbers of the cations and anions in the studied LTILs are 0.77 < t+ < 0.81 and 0.19 < t− < 0.23, respectively. Kowsari et al.7 showed that for MILs 0.49 < t+ < 0.65, and Yeganegi et al.31 showed that for DILs 0.51 < t+ < 0.73. Also, Wu et al.69 reported a t+ value of 0.637 for DIL 3,3-(2,2-oxybis(ethane-2,1diyl))bis(1-allyl-imidazolium)bis(trifluoromethylsulfonyl)amide ([AIOIA][TFSI]), which is in the above range for DILs. The results show that similar to MILs and DILs, cations have a major contribution toward electric current, but this contribution increases from MILs to LTILs. The larger transport numbers of the cations for these LTILs can be related to the larger diffusion coefficients of the cations and also to the larger charges of the cations with respect to those of the anions. To investigate the temperature dependence of the simulated diffusion, viscosity, and electrical conductivity of the studied LTILs, there is a need to perform sufficiently long simulations to achieve an accurate evaluation of these parameters at each temperature. Figure 6 shows the results of such simulations for LTIL-1. The symbols are simulated results and lines are due to fittings with the Vogel−Fulcher−Tammann (VFT) equation.70 This equation is widely used to show the temperature dependence of the transport properties of ILs at T > Tx (where Tx is fragile-to-strong temperature). The VFT equation is given via eq 10 A = A′e B ′ /(T − T0)

(10)

In the VFT equation (eq 10) A can stand for diffusion, viscosity, or electrical conductivity of the IL at a given temperature, T, and the fit parameters are the fragility (B′), which is related to the free activation energy of a fluid, Vogel temperature (T0), and limiting high-temperature quantity (A′). As shown in Figure 6, the VFT equation can fit the simulated transport properties of LTIL-1 in the temperature range of 303 to 453 K well. 3.3. Structural Properties. 3.3.1. Polar and Nonpolar Regions. To visualize the polar−nonpolar nanosegregation within the studied LTIL systems, we have used a color code to show the polar and nonpolar regions. Yellow represents the high-charge-density region, including all atoms belonging to the anions plus all atoms belonging or directly attached to the imidazolium ring. All other atoms (the remainder of the alkyl side chains of the cations), which belong to the low-chargedensity area, are shown in blue and are considered as nonpolar domains. Examples of simulation boxes containing the studied LTIL systems are shown in Figure 7 using the abovementioned color scheme. As shown in this figure, on increasing the length of the alkyl chain on the cation of the LTIL, the nonpolar regions percolate the entire simulation box. This led to the idea of ionic liquids as nanosegregated fluids.70 Lopes and Padua71 used a fully atomistic model to examine a series of dialkylimidazolium hexafluorophosphate ILs with varying alkyl chain lengths. They observed that as the alkyl tails increased in length, the nanodomains of polar and nonpolar regions formed more explicitly in the liquid. It seems that for sufficiently long side chains the nonpolar alkyl chains tend to aggregate, which is in agreement with the many experimental observations.47 3.3.2. Radial Distribution Functions (RDFs) and Spatial Distribution Functions (SDFs). The liquid structure in ILs,

Figure 6. Temperature dependence of the simulated diffusion, viscosity, and electrical conductivity of LTIL-1. The symbols are simulated results and lines are fittings with the VFT equation.

which is as a result of the balance between the long-range electrostatic forces and local van der Waals interactions between the cations and anions, can be investigated using site−site RDFs (g(r)) and also SDFs of anions around the cations. The calculated RDFs for cation−cation (CY−CY), cation−cation (CA−CA), anion−anion, and ring−anion correlations for the studied LTILs at 296 K and 1 atm have 1885


Article


been shown in Figure 8. In Figure 8a, which shows the interactions among the side rings of the cations (CY−CY), there is no a significance difference in the RDFs with increasing length of the alkyl chain. There are distinct oscillations readily visible in these RDFs, indicating long-range spatial correlations between ions. The RDFs of cation−cation (CA−CA), which is the interaction between the middle rings of the cations of LTILs, show broad peaks, which decay slowly (Figure 8b). As shown in this figure, the RDFs for LTIL-1, with three −CH2 units between the rings, show a broad first peak at about 8 Å, but for larger cations, the calculated RDFs lack a well-defined first peak. This is expected because the TILs have bulky cations, which indicates the lack of an organized structure of cations around each other. The less sharp and more complex features of these peaks compared to those of the anion−anion peaks and also the cation−anion RDFs may be related to the accessibility and arrangement of the cations in these ILs, which affect the interactions among the middle rings of the cations. As shown in Figure 8c, anion−anion RDFs show a first peak at around 7.5 Å, with a second peak at distances larger than 15 Å. These distinct and split peaks are a result of sequential ordering induced by the cation−anion pairs. According to this figure, lengthening the alkyl chains between the rings in the cation has no significant effect on the anion−anion RDFs. Figure 8d shows the RDFs of the middle rings of the cations and anions. This figure shows that the distribution of the anions around the middle ring of the cation is insensitive to cation type. This can be attributed to the fact that there is a higher probability of finding anions around the ring planes than alkyl chains. In addition, the existence of the second and third weaker peaks in these RDFs indicate that the structural correlation is longranged and anions and cations distribute homogeneously because of the strong electrostatic interactions. Figure 9 shows the RDFs of anion [NTf2]− around alkyl chain atoms C1, C3, C5, C7, and C9 (for labels refer to those in Figure 1) for LTIL-3 at 296 K and 1 atm. As shown in this figure, on moving from the rings toward the center of the alkyl chain, the peaks become shorter and broader. Again, this result shows that 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. RDFs of anion [NTf2]− around the middle ring hydrogen atoms, that is, H1 and H2, side ring hydrogen atoms, that is, H3, H4, and H5, and middle and side ring hydrogen atoms, that is, H1 and H3, for LTIL-3 at 296 K and 1 atm have been shown in Figure 10. The labels are defined according to those in Figure 1. The distinct and sharp peaks at about 3 and 5 Å in Figure 10a− c indicate strong interactions between the hydrogen atoms of the rings of cations and the fluorine atoms of anions. The existence of these sharp peaks at small distances illustrates strong interactions between these charged sites. Also, the existence of broad peaks even beyond 20 Å is due to the longrange Coulombic interactions. Figure 10a,b shows that the interactions between the anions and hydrogens above the ring planes are stronger than the interactions below the ring planes for both middle and side rings. According to Figure 10b, the interaction between the anion and hydrogen H4, which is toward the outside of the cation, is a little stronger than the interaction with H5, which is toward the inside of the cation. Similar results were observed and reported for imidazolium MILs41 and DILs.31 Also, Figure 10c shows that the interaction between the anion and hydrogen of the middle ring of the

Figure 7. Snapshots of the simulation boxes containing 125 LTILs, using a color code to identify the polar (depicted in yellow) and nonpolar (depicted in blue) domains that are formed in the studied LTILs. 1886


Article


Figure 8. Calculated RDFs for (a) cation−cation (CY−CY), (b) cation−cation (CA−CA), (c) anion−anion, and (d) ring−anion correlations for the studied LTILs at 296 K and 1 atm.

cation−anion RDF exhibits a strong first peak at about 3 Å, indicating a preferential cation−anion interaction over anion− anion and cation−cation interactions, which exhibit first peaks at about 8 or 9 Å. In other words, the first peaks of the cation− cation and anion−anion RDFs are positioned at larger separations. This can be due to the fact that the ions of a given type prefer to be surrounded by ions of opposite charge, to preserve the conditions of local electroneutrality. All RDFs exhibit a second less-intense but broader peak and also a weak third peak at longer distances close to or even beyond 20 Å, which is due to long-range Coulombic interactions. The cation−cation RDFs exhibit a broad first peak at about 9 Å, with a low intensity compared to those of the other two pair types. Alkyl chain−alkyl chain RDFs of the studied LTILs have been shown in Figure 12 at 296 K and 1 atm. The alkyl chain− alkyl chain RDFs of the studied LTILs have a sharp peak at about 5 Å, but there is no obvious second peak, which indicates heterogeneous distributions and aggregation of the alkyl chains. It seems that in the studied LTILs unlike the Coulombic interactions between the rings and with the anions, which are determined mainly by the partial charges on these sites, the interactions between the alkyl chains are mainly due to the short-range attractive van der Waals interactions. In fact, these two kinds of interactions compete with each other under a

Figure 9. RDFs of anion [NTf2]− around alkyl chain atoms C1, C3, C5, C7, and C9 (for labels refer to those in Figure 1) for LTIL-3 at 296 K and 1 atm.

cation is a little stronger than that with the hydrogen of the side rings, but this difference is not significant. Figure 11 shows the center of mass RDFs of the studied LTILs at 296 K and 1 atm. For each of the studied LTILs, the 1887


Article


Figure 10. RDFs of anion [NTf2]− around the (a) middle ring hydrogen atoms, that is, H1 and H2, (b) side ring hydrogen atoms, that is, H3, H4, and H5, and (c) middle and side ring hydrogen atoms, that is, H1 and H3, for LTIL-3 at 296 K and 1 atm. The labels are defined according to Figure 1.

Figure 11. Center of mass RDFs of (a) LTIL-1, (b) LTIL-2, and (c) LTIL-3 at 296 K and 1 atm.

aggregation of the alkyl chain groups. In the other words, with increasing length of the alkyl chain between the rings, the alkyl chain groups have more freedom to be closer to each other, although this aggregation is restricted by geometrical constraint. Three-dimensional SDFs of anion [NTf2]− around cations are calculated using the travis package,72 visualized by the VMD program,73 and shown in Figure 13. In this figure, mauve stands for [NTf2]−, cyan for carbon, blue for nitrogen, and white for hydrogen. As expected, similar to the results obtained for

geometrical constraint between the rings and alkyl chain groups. As expected, the intensity of the first peak increases with increasing length of the alkyl chain between the rings. From LTIL-1 to LTIL-3, the −CH2 groups in the alkyl chains between the rings are comparatively far from the rings and the Coulombic interactions on the rings have less effect on the 1888


Article


Figure 12. Alkyl chain−alkyl chain RDFs of the studied LTILs at 296 K and 1 atm.

MILs41 and DILs,31 the probability of finding anions around the imidazolium rings is higher than that around the alkyl chain and anions are clustered around the rings. For better presentation and clarity, the clustering of anions about the rings has been presented just for the middle rings in Figure 13. As is recognizable in the SDFs, the bulky trications bend to maximize the electrostatic interactions of three imidazolium rings with the anions and also to minimize the interactions between ions with the same charges.

4. CONCLUSIONS MD simulations of three imidazolium-based LTILs were carried out using a nonpolarizable all-atom force field, which is a refined version of the CL&P force field.43 The effects of temperature and alkyl chain length on the densities, MSDs, selfdiffusivities, viscosities, electrical conductivities, transference numbers, and hydrodynamic radii of the cations and anions of the studied LTILs were investigated. The calculated densities of the LTILs are greater than those of both MILs and DILs and decrease with increasing temperature and alkyl chain length between the imidazolium rings. Detailed microscopic structures of the studied LTILs were discussed in terms of RDFs and SDFs. The results show that, similar to that in MILs and DILs, the anions are very well organized around the cationic rings and tend to spend most of their time around the ring compared with alkyl chain atoms. Diffusion coefficients of the studied LTILs are smaller than those of both MILs and DILs, with comparable viscosities. It can be observed that, unlike those of MILs and DILs, the diffusion coefficients of the cations and anions of the studied LTILs increase with increasing length of the alkyl chain between the rings from LTIL-1 to LTIL-2 but decrease for LTIL-3, which is in a good agreement with the trend of viscosity data.13 The calculated MSDs and diffusion coefficients of the alkyl chain, middle ring, and side ring in each of the studied LTILs are similar, which indicates that the motion of these groups is highly correlated to each other. The trend of the molar electrical conductivity of the studied LTILs is in contrast to the viscosity values of these ILs. LTIL-2, which has the lowest viscosity, has a greater molar electrical conductivity. The large transference numbers of the cations (0.77 < t+ < 0.81) show that, similar to that in MILs and DILs,

Figure 13. SDFs of anion [NTf2]− around the cations: mauve, [NTf2]−; cyan, carbon; blue, nitrogen; white, hydrogen.

cations have a major contribution to electric current, and this contribution is more than that in both MILs and DILs.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b10766. Electrostatic charges for the atoms of the studied LTILs obtained at the B3LYP/6-31++G(d,p) level of theory by the CHELPG method; evolution of parameter β with time for the middle ring, side ring, alkyl chain, and anion 1889


Article


■

(14) Danna, F.; Ferrante, F.; Noto, R. Geminal Ionic Liquids: A Combined Approach to Investigate Their Three-Dimensional Organization. Chem. Eur. J. 2009, 15, 13059−13068. (15) Engel, R.; Cohen, J. I. Organic Polycationic Salts-Syntheses and Applications. Curr. Org. Chem. 2002, 6, 1453−1467. (16) Wishart, J. F.; Lall-Ramnarine, S. I.; Raju, R.; Scumpia, A.; Bellevue, S.; Ragbir, R.; Engel, R. Effects of Functional Group Substitution on Electron Spectra and Solvation Dynamics in a Family of Ionic Liquids. Radiat. Phys. Chem. 2005, 72, 99−104. (17) Qi, M.; Armstrong, D. W. Dicationic Ionic Liquid Stationary Phase for GC-MS Analysis of Volatile Compounds in Herbal Plants. Anal. Bioanal. Chem. 2007, 388, 889−899. (18) Han, X.; Armstrong, D. W. Using Geminal Dicationic Ionic Liquids as Solvents for High-Temperature Organic Reactions. Org. Lett. 2005, 7, 4205−4208. (19) Jin, C. M.; Ye, C.; Phillips, B. S.; Zabinski, J. S.; Liu, X.; Liu, W.; Shreeve, J. M. Glycol Functionalized Dicationic Ionic Liquids with Alkyl or Polyfluoro Alkyl Substituents as High Temperature Lubricants. J. Mater. Chem. 2006, 16, 1529−1535. (20) Zhang, Z.; Zhoua, H.; Yanga, L.; Tachibana, K.; Kamijima, K.; Xu. Asymmetrical Dicationic Ionic Liquids Based on both Imidazolium and Aliphatic Ammonium as Potential Electrolyte Additives Applied to Lithium Secondary Batteries. Electrochim. Acta 2008, 53, 4833−4838. (21) Kim, J. Y.; Kim, T. H.; Kim, D. Y.; Park, N. -G; Ahn, K. D. Novel Thixotropic Gel Electrolytes Based on Dicationic Bisimidazolium Salts for Quasi-Solid-State Dye-Sensitized Solar Cells. J. Power Sources 2008, 175, 692−697. (22) Han, X.; Armstrong, D. W. Ionic Liquids in Separations. Acc. Chem. Res. 2007, 40, 1079−1086. (23) Mutelet, F.; Moise, J. C.; Skrzypczak, A. Evaluation of the Performance of Trigeminal Tricationic Ionic Liquids for Separation Problems. J. Chem. Eng. Data 2012, 57, 918−927. (24) Payagala, T.; Zhang, Y.; Wanigasekara, E.; Huang, K.; Breitbach, Z. S.; Sharma, P. S.; Sidisky, L. M.; Armstrong, D. W. Trigonal Tricationic Ionic Liquids: A Generation of Gas Chromatographic Stationary Phases. Anal. Chem. 2009, 81, 160−173. (25) Breitbach, Z. S.; Warnke, M. M.; Wanigasekara, E.; Zhang, X.; Armstrong, D. W. Evaluation of Flexible Linear Tricationic Salts as Gas-Phase Ion-Pairing Reagents for the Detection of Divalent Anions in Positive Mode ESI-MS. Anal. Chem. 2008, 80, 8828−8834. (26) Bodo, E.; Caminiti, R. The Structure of Geminal Imidazolium Bis(trifluoromethylsulfonyl) amide Ionic Liquids: A Theoretical Study of the Gas Phase Ionic Complexes. J. Phys. Chem. A 2010, 114, 12506−12512. (27) Bodo, E.; Chiricotto, M.; Caminiti, R. Structure of Geminal Imidazolium Bis(trifluoromethylsulfonyl)imide Dicationic Ionic Liquids: A Theoretical Study of the Liquid Phase. J. Phys. Chem. B 2011, 115, 14341−14347. (28) Vitorino, J.; Leal, P.; Licence, P.; Lovelock, K. R. J.; Gooden, P. N.; Piedade, M. E. M.; Shimizu, K.; Rebelo, L. P. N.; Lopes, J. N. C. Vaporization of a Dicationic Ionic Liquid Revisited. ChemPhysChem 2010, 11, 3673−3677. (29) Bhargava, B. L.; Klein, M. L. Formation of Interconnected Aggregates in Aqueous Dicationic Ionic Liquid Solutions. J. Chem. Theory Comput. 2010, 6, 873−879. (30) Bhargava, B. L.; Klein, M. L. Nanoscale Organization in Aqueous Dicationic Ionic Liquid Solutions. J. Phys. Chem. B 2011, 115, 10439−10446. (31) Yeganegi, S.; Soltanabadi, A.; Farmanzadeh, D. Molecular Dynamic Simulation of Dicationic Ionic Liquids: Effects of Anions and Alkyl Chain Length on Liquid Structure and Diffusion. J. Phys. Chem. B 2012, 116, 11517−11526. (32) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (33) Lopes, J. N. C.; Deschamps, J.; Padua, A. A. H. Modeling Ionic Liquids using a Systematic All-Atom Force Field. J. Phys. Chem. B 2004, 108, 2038−3047.

for LTIL-1, LTIL-2, and LTIL-3 at 303 K and 1 atm (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: +98-313-7934942. Fax: +98-313-668-9732. ORCID

Majid Moosavi: 0000-0003-3136-3464 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported financially by the Research Council of the University of Isfahan. The authors are indebted to the Institute for Studies in Theoretical Physics and Mathematics (IPM) for providing the HPC Cluster to perform a part of the simulations. Also, they thank Prof. H. Sabzyan, M. Bahrami, and T. Sedghamiz for their guidance.

■

REFERENCES

(1) Earle, M. J.; Esperança, J. M. S. S.; Gilea, M. A.; Canongia Lopes, J. N.; Rebelo, L. P. N.; Magee, J. W.; Seddon, K. R.; Widegren, J. A. The Distillation and Volatility of Ionic Liquids. Nature 2006, 439, 831−834. (2) MacFarlane, D. R.; Forsyth, M.; Izgorodina, E. I.; Abbott, A. P.; Annat, G.; Fraser, K. On the Concept of Ionicity in Ionic Liquids. Phys. Chem. Chem. Phys. 2009, 11, 4962−4967. (3) Weingarth, D.; Foelske-Schmitz, A.; Koetz, R. PTFE Bound Activated Carbon − A Quasi Reference Electrode for Ionic Liquids and Its Application. J. Power Sources 2013, 225, 84−88. (4) Wassercheid, P.; Welton, T. Ionic Liquids in Synthesis, 2nd ed.; Wiley-VCH: New York, 2008. (5) Welton, T. Room-Temperature Ionic Liquids. Solvents for Synthesis and Catalysis. Chem. Rev. 1999, 99, 2071−2083. (6) Tokuda, H.; Hayamizu, K.; Ishii, K.; Susan, M. A. B. H.; Watanabe, M. Physicochemical Properties and Structures of Room Temperature Ionic Liquids. 2. Variation of Alkyl Chain Length in Imidazolium Cation. J. Phys. Chem. B 2005, 109, 6103−6110. (7) Kowsari, M. H.; Alavi, S.; Ashrafizaadeh, M.; Najafi, B. Molecular Dynamics Simulation of Imidazolium-Based Ionic Liquids. I. Dynamics and Diffusion Coefficient. J. Chem. Phys. 2008, 129, No. 224508. (8) Paredes, X.; Fernández, J.; Pádua, A. H.; Malfreyt, P.; Malberg, F.; Kirchner, B.; Pensado, A. S. Bulk and Liquid−Vapor Interface of Pyrrolidinium-Based Ionic Liquids: A Molecular Simulation Study. J. Phys. Chem. B 2014, 118, 731−742. (9) Sharma, P. S.; Payagala, T.; Wanigasekara, E.; Wijeratne, A. B.; Huang, J.; Armstrong, D. W. Trigonal Tricationic Ionic Liquids: Molecular Engineering of Trications to Control Physicochemical Properties. Chem. Mater. 2008, 20, 4182−4184. (10) Shirota, H.; Mandai, T.; Fukazawa, H.; Kato, T. Comparison between Dicationic and Monocationic Ionic Liquids: Liquid Density, Thermal Properties, Surface Tension, and Shear Viscosity. J. Chem. Eng. Data 2011, 56, 2453−2459. (11) Anderson, J. L.; Ding, R.; Ellern, A.; Armstrong, D. W. Structure and Properties of High Stability Geminal Dicationic Ionic Liquids. J. Am. Chem. Soc. 2005, 127, 593−604. (12) Payagala, T.; Huang, J.; Breitbach, Z. S.; Sharma, P. S.; Armstrong, D. W. Unsymmetrical Dicationic Ionic Liquids: Manipulation of Physicochemical Properties using Specific Structural Architectures. Chem. Mater. 2007, 19, 5848−5850. (13) Wanigasekara, E.; Zhang, X.; Nanayakkara, Y.; Payagala, T.; Moon, H.; Armstrong, D. W. Linear Tricationic Room-Temperature Ionic Liquids: Synthesis, Physiochemical Properties, and Electrowetting Properties. ACS Appl. Mater. Interfaces 2009, 1, 2126−2133. 1890


Article


(55) Yeh, I.-C.; Hummer, G. System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions. J. Phys. Chem. B 2004, 108, 15873−15879. (56) Prado, C. E. R.; Freitas, L. C. G. Molecular Dynamics Simulation of the Room-Temperature Ionic Liquid 1-Butyl-3Methylimidazolium Tetrafluoroborate. J. Mol. Struc. THEOCHEM 2007, 84, 93−100. (57) Zhong, X.; Liu, Z.; Cao, D. Improved Classical United-Atom Force Field for Imidazolium-Based Ionic Liquids: Tetrafluoroborate, Hexafluorophosphate, Methylsulfate, Trifluoromethylsulfonate, Acetate, Trifluoroacetate, and Bis(trifluoromethylsulfonyl)amide. J. Phys. Chem. B 2011, 115, 10027−10040. (58) Alavi, S.; Thompson, D. L. Simulations of the Solid, Liquid, and Melting of 1-n-Butyl-4-amino-1,2,4-triazolium Bromide. J. Phys. Chem. B 2005, 109, 18127−18134. (59) Chaban, V. V.; Voroshylovab, I. V.; Kalugin, O. N. A New Force Field Model for the Simulation of Transport Properties of Imidazolium-Based Ionic Liquids. Phys. Chem. Chem. Phys. 2011, 13, 7910−7920. (60) Shimizu, K.; Tariq, M.; Gomes, M. F. C.; Rebelo, L. P. N.; Lopes, J. N. C. Assessing the Dispersive and Electrostatic Components of the Cohesive Energy of Ionic Liquids Using Molecular Dynamics Simulations and Molar Refraction Data. J. Phys. Chem. B 2010, 114, 5831−5834. (61) Tariq, M.; Serro, A. P.; Mata, J. L.; Saramago, B.; Esperança, J. M.; Lopes, J. N. C.; Rebelo, L. P. N. High-Temperature Surface Tension and Density Measurements of 1-alkyl-3-Methylimidazolium Bistriflamide Ionic Liquids. Fluid Phase Equilib. 2010, 294, 131−138. (62) Moosavi, M.; Khashei, F.; Sharifi, A.; Mirzaei, M. The Effects of Temperature and Alkyl Chain Length on the Density and Surface Tension of the Imidazolium-Based Geminal Dicationic Ionic Liquids. J. Chem. Thermodyn. 2017, 107, 1−7. (63) Hasse, B.; Lehmann, J.; Assenbaum, D.; Wasserscheid, P.; Leipertz, A.; Fröba, A. P. Viscosity, Interfacial Tension, Density, and Refractive Index of Ionic Liquids [EMIM][MeSO3], [EMIM][MeOHPO2], [EMIM][OcSO4], and [BBIM][NTf2] in Dependence on Temperature at Atmospheric Pressure. J. Chem. Eng. Data 2009, 54, 2576−2583. (64) Jacquemin, J.; Husson, P.; Paduaa, A. A. H.; Majera, V. Density and Viscosity of Several Pure and Water-saturated Ionic Liquids. Green Chem. 2006, 8, 172−180. (65) Kowsari, M. H.; Alavi, S.; Ashrafizaadeh, M.; Najafi, B. Molecular Dynamics Simulation of Imidazolium-Based Ionic Liquids. II. Transport Coefficients. J. Chem. Phys. 2009, 130, No. 014703. (66) Van-Oanh, N. T.; Houriez, C.; Rousseau, B. Viscosity of the 1Ethyl-3-Methylimidazolium Bis(Trifluoromethylsulfonyl)imide Ionic Liquid from Equilibrium and Nonequilibrium Molecular Dynamics. Phys. Chem. Chem. Phys. 2010, 12, 930−936. (67) Jiang, W.; Yan, T.; Wang, Y.; Voth, G. A. Molecular Dynamics Simulation of the Energetic Room-Temperature Ionic Liquid, 1Hydroxyethyl-4-amino-1,2,4-triazolium Nitrate (HEATN). J. Phys. Chem. B 2008, 112, 3121−3131. (68) Atkins, P. W.; Julio De, P. Physical Chemistry, 9th ed.; Oxford University Press: U.K., 2010. (69) Wu, T. Y.; Su, S. G.; Wang, Y. S.; Lin, Y. C.; Chang, J. K.; Kuo, C. W.; Sun, I. W. Diffusion Coefficients, Spin-Lattice Relaxation Times, and Chemical Shift Variations of NMR Spectra in LiTFSIDoped Ether- and Allyl-Functionalized Dicationic Ionic Liquids. J. Taiwan Inst. Chem. Eng. 2016, 60, 138−150. (70) Shimizu, K.; Tariq, M.; Rebelo, L. P. N.; Lopes, J. N. C. Binary Mixtures of Ionic Liquids with a Common Ion Revisited: A Molecular Dynamics Simulation Study. J. Mol. Liq. 2010, 153, 52−56. (71) Lopes, J. N. A. C.; Padua, A. A. H. Nanostructural Organization in Ionic Liquids. J. Phys. Chem. B 2006, 110, 3330−3335. (72) Brehm, M.; Kirchner, B. TRAVIS - A Free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories. J. Chem. Inf. Model 2011, 51, 2007−2023.

(34) Lopes, J. N. C.; Padua, A. A. H. Molecular Force Field for Ionic Liquids Composed of Triflate or Bistriflylimide Anions. J. Phys. Chem. B 2004, 108, 16893−16898. (35) Lopes, J. N. C.; Padua, A. A. H. Molecular Force Field for Ionic Liquids III: Imidazolium, Pyridinium, and Phosphonium Cations; Chloride, Bromide, and Dicyanamide Anions. J. Phys. Chem. B 2006, 110, 19586−19592. (36) Lopes, J. N. C.; Padua, A. A. H.; Shimizu, K. Molecular Force Field for Ionic Liquids IV: Trialkylimidazolium and Alkoxycarbonylimidazolium Cations; Alkylsulfonate and Alkylsulfate Anions. J. Phys. Chem. B 2008, 112, 5039−5046. (37) Chaban, V. Polarizability versus Mobility: Atomistic Force Field for Ionic Liquids. Phys. Chem. Chem. Phys. 2011, 13, 16055−16062. (38) Salanne, M. Simulations of Room Temperature Ionic Liquids: from Polarizable to Coarse-Grained Force Fields. Phys. Chem. Chem. Phys. 2015, 17, 14270−14279. (39) Leontyev, I.; Stuchebrukhov, A. Accounting for Electronic Polarization in Non-Polarizable Force Fields. Phys. Chem. Chem. Phys. 2011, 13, 2613−2626. (40) Schmidt, J.; Krekeler, C.; Dommert, F.; Zhao, Y.; Berger, R.; Site, L. D.; Holm, C. Ionic Charge Reduction and Atomic Partial Charges from First-Principles Calculations of 1,3-Dimethylimidazolium Chloride. J. Phys. Chem. B 2010, 114, 6150−6155. (41) Ghatee, M. H.; Zolghadr, A. R.; Moosavi, F.; Ansari, Y. Studies of Structural, Dynamical, and Interfacial Properties of 1-alkyl-3methylimidazolium Iodide Ionic Liquids by Molecular Dynamics Simulation. J. Chem. Phys. 2012, 136, 124706−14. (42) Bhargava, B. L.; Klein, M. L. Refined Potential Model for Atomistic Simulations of Ionic Liquid [bmim][PF6]. J. Chem. Phys. 2007, 127, No. 114510. (43) Tsuzuki, S.; Shinoda, W.; Saito, H.; Mikami, M.; Tokuda, H.; Watanabe, M. Molecular Dynamics Simulations of Ionic Liquids: Cation and Anion Dependence of Self-Diffusion Coefficients of Ions. J. Phys. Chem. B 2009, 113, 10641−10649. (44) Smith, W.; Todorov, I. DL_POLY Classic Molecular Simulation Package; Daresbury Laboratory: Cheshire, U.K., 2012. (45) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−525. (46) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (47) Gardas, R. L.; Coutinho, J. A. P. Group Contribution Methods for the Prediction of Thermophysical and Transport Properties of Ionic Liquids. AIChE J. 2009, 55, 1274−1290. (48) Coutinho, J. A. P.; Carvalho, P. J.; Oliveira, N. M. C. Predictive Methods for the Estimation of Thermophysical Properties of Ionic Liquids. RSC Adv. 2012, 2, 7322−7346. (49) Paduszyński, K.; Domanska, U. A New Group Contribution Method for Prediction of Density of Pure Ionic Liquids over a Wide Range of Temperature and Pressure. Ind. Eng. Chem. Res. 2012, 51, 591−604. (50) Yeganegi, S.; Sokhanvarana, V.; Soltanabadia, A. Study of Thermodynamic Properties of Imidazolium Based Ionic Liquids and Investigation of the Alkyl Chain Length Effect by Molecular Dynamics Simulation. Mol. Simul. 2013, 39, 1070−1078. (51) Ghatee, M. H.; Moosavi, F.; Zolghadr, A. R.; Jahromi, R. Critical-Point Temperature of Ionic Liquids from Surface Tension at Liquid−Vapor Equilibrium and the Correlation with the Interaction Energy. Ind. Eng. Chem. Res. 2010, 49, 12696−12701. (52) Harris, K. R.; Woolf, L. A.; Kanakubo, M. Temperature and Pressure Dependence of the Viscosity of the Ionic Liquid 1-Butyl-3methylimidazolium Hexafluorophosphate. J. Chem. Eng. Data 2005, 50, 1777−1782. (53) Trivedi, S.; Pandey, S. Densities of 1-Butyl-3-methylimidazolium Hexafluorophosphate + Poly (Ethylene Glycol) in the Temperature Range (283.15 to 363.15) K. J. Chem. Eng. Data 2011, 56, 2168−2174. (54) Chaban, V. V.; Prezhdo, O. V. Computationally Efficient Prediction of Ionic Liquid Properties. J. Phys. Chem. Lett. 2014, 5, 1973−1977. 1891


Article

The Journal of Physical Chemistry B (73) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38.

1892


Tricationic Ionic Liquids: Structural and Dynamical Properties via

Recommend Documents