Molecular Modeling of Imidazolium-Based [Tf2Nâ] Ionic Liquids

J. Phys. Chem. B 2009, 113, 7211–7224

7211

Molecular Modeling of Imidazolium-Based [Tf2N-] Ionic Liquids: Microscopic Structure, Thermodynamic and Dynamic Properties, and Segmental Dynamics Georgia-Evangelia Logotheti,†,§ Javier Ramos,‡ and Ioannis G. Economou*,‡ Molecular Thermodynamics and Modeling of Materials Laboratory, Institute of Physical Chemistry, National Center for Scientific Research “Demokritos”, GR-153 10 Aghia ParaskeVi Attikis, Greece, and Department of Macromolecular Physics, Instituto de Estructura de la Materia - CSIC, Serrano 113bis, ES-28006, Madrid, Spain ReceiVed: August 5, 2008; ReVised Manuscript ReceiVed: February 24, 2009

The microscopic structure, thermodynamic properties, local segmental dynamics, and self-diffusion coefficients of three ionic liquids (ILs) with a common anion, namely, the bis(trifluoromethylsulfonyl) imide ([Tf2N-]), and imidazolium-based cations that differ in the alkyl tail length, namely, the 1-butyl-3-methylimidazolium ([C4mim+]), the 1-hexyl-3-methylimidazolium ([C6mim+]), and the 1-octyl-3-methylimidazolium ([C8mim+]), are calculated over the temperature range of 298.15-333.15 K and pressure range of 0.1-60 MPa. Quantum calculations based on density functional theory are performed on isolated ion pairs, and minimum energy conformers are identified. Electronic density results are used to estimate the electrostatic potential of a molecular force field that is used subsequently for long molecular dynamics (MD) simulations of bulk ILs. Thermodynamic properties calculated from MD are shown to be in excellent agreement for the bulk density and good agreement for derivative properties when compared to experimental data. The new force field is an improvement over earlier ones for the same ILs. The microscopic structure as expressed through the radial distribution function is thoroughly calculated, and it is shown that the bulk structure characteristics are very similar to those obtained from the quantum calculations on isolated ion pairs. The segmental dynamics expressed in terms of bond and torsion angle decorrelation is shown to assume a broad range of characteristic times. Molecular segments in the alkyl tail of the cations are significantly faster than segments in the vicinity of the imidazolium ring. Finally, the new force field predicts accurately the self-diffusion coefficients of the cations and the anions over the entire temperature range examined, thus confirming its validity for a broad range of physical properties. 1. Introduction Ionic Liquids (ILs) have received tremendous attention by both the academic and industrial community in recent years because of the emphasis on the development of green-chemistrybased processes and of a number of unique physicochemical properties that make them appropriate candidates for various scientific and technological applications, such as reaction media, catalysis, separation processes, and even gas storage.1-5 These substances usually are made of an asymmetric organic cation (often a derivative of pyridinium, imidazolium, alkylphosphonium, or alkylammonium ion) and an inorganic anion, such as tetrafluoroborate (BF4-), hexafluorophosphate (PF6-), chloride (Cl-), or bis(trifluoromethylsulfonyl) imide (Tf2N-), among others. ILs are nonvolatile salts that are liquids over a broad temperature range. They are good solvents for a number of organic and even inorganic substances and are not flammable and are odorless. Although a lot of effort has been devoted to the experimental determination of the physical properties of various ILs, the available data for many properties are still limited and, in some cases, questionable or conflicting.6 In this respect, accurate theoretical and/or computational methods are * To whom correspondence should be addressed. E-mail: economou@ chem.demokritos.gr. Tel.: +30 210 6503963. † National Center for Scientific Research “Demokritos”. ‡ Instituto de Estructura de la Materia - CSIC. § Current address: Hyperion Systems Engineering, Modeling and Simulation Solutions, 10-12 Doryleou Street, GR-11521 Athens, Greece.

very valuable since they can provide an independent approach for the estimation of properties of interest. Molecular simulation using a realistic representation of interand intramolecular interactions (force fields) is a powerful tool for the elucidation of molecular mechanisms and accurate prediction of microscopic structure and macroscopic properties of ILs.7,8 The accurate parametrization of the force field is based on experimental measurements (i.e., spectroscopic studies, calorimetric data, etc.) and/or quantum mechanics calculations. A number of papers have appeared recently concerning the structure, thermodynamics, and dynamics of imidazolium-based ILs with various anions as elucidated from molecular dynamics (MD) simulations.9-14 IUPAC has selected [Tf2N-] as a benchmark anion in ILs for extensive investigation. As a result, [Tf2N-] ILs have received significant attention in terms of experimental investigation15-17 and molecular modeling17-27 of their structure, thermodynamic (density, enthalpy of vaporization, etc.), and transport (self-diffusion coefficient, viscosity, etc.) properties. Canongia Lopes and Pa´dua19 parametrized a molecular force field for [Tf2N-] in the framework of OPLS-AA/AMBER force fields using quantum mechanics calculations and validated it against experimental structure and thermodynamic properties.19,20 Recently, Ludwig and co-workers reparameterized further this force field to improve agreement with experimental thermodynamic and dynamic data for [Tf2N-] ILs.21,22 Besides their unique bulk thermodynamic properties compared to other liquids, ILs exhibit also unique heterogeneous

10.1021/jp806999y CCC: $40.75  2009 American Chemical Society Published on Web 04/30/2009

7212

J. Phys. Chem. B, Vol. 113, No. 20, 2009

Logotheti et al.

Figure 1. Atom numbering for the ILs examined in this work.

behavior at the microscopic level, as elucidated recently by detailed atomistic23 and more coarse-grained24,25 molecular simulation and confirmed by X-ray diffraction experiments.28 On the basis of the majority of the experimental and simulation data discussed above, it is clear that the sole variability of the ion identity (anion or cation) in a given IL may lead to remarkable changes in their physicochemical properties. Thus, an understanding at the molecular level is highly desirable, helping lead to a rational design oriented to a specific application. This work focuses on three ILs with a common anion, namely, the [Tf2N-]. The three ILs that are examined here include 1-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl) imide ([C4mim+][Tf2N-]), 1-hexyl3-methylimidazolium bis(trifluoromethylsulfonyl) imide ([C6mim+][Tf2N-]), and 1-octyl-3-methylimidazolium bis(trifluoromethylsulfonyl) imide ([C8mim+][Tf2N-]). In this way, it is possible to investigate thoroughly the structural and dynamic properties of ILs as a function of the alkyl chain length at a range of temperatures and pressures. Quantum calculations based on density functional theory (DFT) are performed initially on isolated ion pairs. Minimum-energy conformers are identified, and results are used to parametrize a molecular force field that is used subsequently for bulk MD simulations at various temperature and pressure values. This approach differs from earlier attempts to parametrize molecular force fields for ILs based on quantum calculations for isolated cations and anions, respectively.19 The approach used here is more realistic in the calculation of the electrostatic terms of the force fields. Thermodynamic properties such as density, thermal expansion coefficient, and isothermal compressibility are calculated and compared against literature experimental data. In addition, the microscopic structure of the ILs in terms of the radial distribution function is calculated. Moreover, very long MD simulations on

Figure 2. Electrostatic potential mapped on the electronic density for the [C6mim+] cation and [Tf2N-] anion at the B3LYP/6-311+G* level. The isosurface value is 0.02 e. The electrostatic potential ranges from -330 (red) to -125 kJ/mol (yellow) for [Tf2N-] and from 85 (yellow) to 770 kJ/mol (blue) for [C6mim+].

the order of 50 ns allow reliable calculation of the local (segmental) molecular dynamics as well as dynamic properties such as the self-diffusion coefficient. In all cases, simulation results are in very good agreement with experiments and, in most cases, are more accurate than predictions from earlier simulation work.19,21 2. Simulation Details 2.1. Quantum Calculations. Quantum calculations were performed on isolated [Cnmim+][Tf2N-] ion pairs with n ) 4, 6, or 8, based on DFT theory in order to establish the relative position and orientation of the anion with respect to the cation and the atomic electrostatic charges on the ions to parametrize

Molecular Modeling of Imidazolium-Based [Tf2N-]

J. Phys. Chem. B, Vol. 113, No. 20, 2009 7213

Figure 3. Relative electronic energies of the isolated ion pairs optimized at the B3LYP/6-311+G* level. Energies are given in kJ/mol. Conformer III has been arbitrarily chosen as the reference; thus, relative energies are calculated by ∆Ei ) Ei - E3. Energies with ZPE corrections are given in parentheses.

TABLE 1: Selected Geometrical Parameters for Various [Cnmim+][Tf2N-] Ion Pair Conformers from DFT Calculations conformer

C2-N4 (Å)

C2-H1 (Å)

C2-O1 (Å)

C2-O2 (Å) +

I II III IV V VI

3.117 3.117 3.117 4.687 3.181 3.467

1.089 1.089 1.089 1.083 1.089 1.083

3.636 3.647 3.642 5.366 3.495 3.043

I II III IV V VI

3.125 3.124 3.122 4.667 3.178 3.484

1.087 1.089 1.089 1.083 1.089 1.083

3.577 3.645 3.643 5.309 3.477 3.041

I II III IV V VI

3.127 3.121 3.122 4.694 3.178 3.498

1.088 1.089 1.089 1.083 1.089 1.083

3.630 3.663 3.660 5.393 3.486 3.037

[C4min ][Tf2N ] 5.316 5.324 5.320 2.982 5.282 4.918 [C6min+][Tf2N-] 5.223 5.344 5.332 2.984 5.269 4.944 [C8min+][Tf2N-] 5.296 5.345 5.343 2.983 5.275 4.926

the atomistic force field to be used subsequently in the classical MD simulations. All energies and geometries of the isolated ion pairs were calculated using the B3LYP hybrid DFT method as implemented in Gaussian03.29 This method consists of two different functionals, the Becke three-parameter hybrid exchange functional30 and the Lee-Yang-Parr correlation functional.31 Both potentials use nonlocal exchange-correlation corrections for the density. The electronic configurations of the atoms were described by the 6-311+G* basis set. Atomic charges of different ionic pair conformers were calculated by electrostatic surface potential fits, using the CHelpG procedure,32 to electron densities obtained at the B3LYP/6-311+G* level. In their force field development for [Tf2N-] ILs, Canongia Lopes and Pa´dua19 performed

C2-O3 (Å)

C2-O4 (Å)

X-H1-N4-S2 (deg)

3.364 3.360 3.361 5.480 3.615 3.233

4.724 4.729 4.727 3.007 4.983 3.488

138.4 136.1 137.1 76.4 -0.8 157.0

3.300 3.436 3.403 5.477 3.618 3.250

4.421 4.818 4.774 3.009 4.981 3.492

122.8 164.0 141.6 78.3 14.2 160.4

3.340 3.382 3.385 5.483 3.614 3.184

4.596 4.760 4.758 3.004 4.969 3.525

125.2 138.9 140.3 76.8 11.7 135.9

-

quantum mechanics calculations at the MP2/cc-pVTZ(-f) level using HF/6-31G(d) to optimize the geometries. This method is clearly inferior to the method used here for both optimization and energy calculation. It should be mentioned that electrostatic charges are very dependent on the geometrical structures; thus, in order to have a set of “good” electrostatic charges, it is important to have better geometries than energies. In the literature, most authors provide atomic charges for the isolated cation and for the isolated anion but not for the cation-anion ion pair. In this way, the force field is transferable between different ILs. Here, the atomic charges for the most stable ion pair conformers are calculated and averaged over all conformers, thus providing a new way to calculate the charges

7214


Logotheti et al. TABLE 2: Mass Density of [C4mim+][Tf2N-], [C6mim+][Tf2N-], and [C8mim+][Tf2N-] at Various Temperatures and Pressures from NPT MD Simulations mass density (g/cm3)

T (K) 0.1 MPa

26 MPa +

[C4mim ][Tf2N ] 1.452 ( 0.003 1.467 ( 0.001 1.439 ( 0.003 1.456 ( 0.004 1.428 ( 0.003 1.447 ( 0.003 1.418 ( 0.003 1.437 ( 0.003 [C6mim+][Tf2N-] 1.378 ( 0.005 1.392 ( 0.007 1.369 ( 0.002 1.384 ( 0.005 1.358 ( 0.003 1.377 ( 0.005 1.343 ( 0.002 1.363 ( 0.005 [C8mim+][Tf2N-] 1.331 ( 0.003 1.321 ( 0.002 1.314 ( 0.003 1.304 ( 0.003

298.15 308.15 318.15 328.15 Figure 4. Experimental data (lines) and MD calculations (points) for the mass densities of [C4mim+][TF2N-] (black), [C6mim+][TF2N-] (red), and [C8mim+][TF2N-] (green). Squares indicate calculations at 0.1 MPa, circles those at 26 MPa, and upper triangles those at 60 MPa, respectively.

for the force field. Clearly, this approach results in a less transferable set of parameters (DFT calculations are needed for each conformer and for each IL), but it is expected to be more “realistic” in the determination of electrostatic terms. A clear justification for this is provided by the excellent prediction of thermodynamic and dynamic properties of ILs shown below. In addition, the electrostatic potential mapped on the electronic density was calculated for the isolated [Cnmim+] cations and [Tf2N-] anion in order to find the most positive and negative areas of the two species. Electrostatic surfaces were found to be very similar for all cations. Different initial guess conformations for the three isolated ion pairs were prepared by matching the most positive areas of the cation with the more negative areas in the anion. Then, initial geometries were optimized using the B3LYP/6-311+G* method described above. Vibrational analyses on all calculated structures revealed a lack of imaginary frequencies, confirming the presence of true minima on the potential energy surface. 2.2. Classical MD Simulations. MD simulations were performed on systems consisting of 100 [Cnmim+][Tf2N-] ion pairs with n ) 4, 6, or 8. For the construction of the initial structure, the amorphous builder module implemented in MAPS 3.0 software package was used.33 The initial structure was relaxed using a conjugate gradient energy minimization procedure, also available in MAPS 3.0, until atomic overlaps were considerably reduced. Subsequently, NPT MD simulations were carried out over a temperature range from 298.15 to 333.15 K and a pressure range from 0.1 to 60 MPa using the NAMD simulation package.34 Typical simulation runs were on the order of 20-50 ns, depending on the size of the alkyl chain and the temperature, in order to capture accurately the dynamics of the ILs. A very efficient full electrostatics algorithm based on the Particle Mesh Ewald (PME) method was used to deal with the full electrostatic interactions. The equations of motion were integrated using the velocity Verlet method. To speed up the MD simulation, a multiple time step algorithm was used, namely, the reversible reference system propagator algorithm (rRESPA).35,36 This algorithm decouples slow and fast degrees of freedom using different Liouville operators for each category, thus allowing the use of more than one time step for the time integration of the different degrees of freedom. In this work, the time step for the bonded interactions was 1 fs, while for the short-ranged nonbonded van der Waals interactions and the full electrostatic interactions, larger time steps of 2 and 4 fs, respectively, were used. The temperature was controlled using a Langevin piston method, while the Nose´-Hoover barostat was used to control pressure.37 Potential tails were brought smoothly to zero using a switching function.38

60 MPa

-

298.15 308.15 318.15 333.15 298.15 308.15 318.15 328.15

1.481 ( 0.003 1.480 ( 0.002 1.468 ( 0.002 1.460 ( 0.004 1.412 ( 0.006 1.403 ( 0.005 1.395 ( 0.005 1.382 ( 0.005

A detailed all-atom force field was used to model ILs given by the expression

U ) Ubonded + Unonbonded )

∑ kb(b - bo)2 + ∑ kθ(θ -

bonds

angles

4

θo)2 +

∑ ∑ kx[1 + cos(nx + δ)] + ∑

dihedrals n)1 n-1

ψo)2 +

n

∑∑ i)1

kψ(ψ -

impropers

j>i

{ [( ) ( ) ] 4εij

σij rij

12

-

σij rij

6

+

qiqj 4ε0rij

}

(1)

where b, θ, and ψ denote bond length, bond angle, and improper angle, respectively. A subscript o stands for the equilibrium value, while ε and σ correspond to the Lennard-Jones parameters. Partial charges are denoted by qi, and ε0 stands for the permittivity in vacuum. All bonded and Lennard-Jones parameters were taken from Maginn and co-workers39,40 for the cations and from Canongia Lopes and Pa´dua19 for the anion. The partial charges for the electrostatic interactions were taken from the quantum DFT calculations presented here. All of the force field parameter values that were used in this work for the various atoms are provided as Supporting Information. In Figure 1, a schematic representation of the three ILs is provided together with appropriate atom numbering that is used throughout this paper. The initial configurations were equilibrated for a period ranging from 1 to 5 ns, depending on the simulation conditions, the initial structure, and the length of the cation’s alkyl chain. Stabilization of the running average for the density and the various components of the energy were used as a criterion for equilibration. Configurations generated during the production stage of each simulation were recorded every 1 ps in order to track accurately the motions that dominate the short-time segmental dynamics. Statistical uncertainties in the thermodynamic properties and relaxation times reported in this work were calculated based on block average analysis. Parallel DFT calculations were carried out at the CIEMAT center in Spain. Most MD simulations were run on a Xeon cluster at the CTI-CSIC Centre in Spain. 3. Results and Discussion 3.1. DFT Calculations. 3.1.1. Electrostatic Potential of Isolated Ion Pairs. The electrostatic potential mapped on the electronic density evaluated



TABLE 3: Thermal Expansion Coefficient, aP, of [C4mim+][Tf2N-], [C6mim+][Tf2N-,] and [C8mim+][Tf2N-] at 298.15 K from Experimental Measurements and NPT MD Simulations aP (×104 K-1)

ionic liquid 0.1 MPa expt +

-

[C4mim ][Tf2N ] [C6mim+][Tf2N-] [C8mim+][Tf2N-]

44

26 MPa MD

6.36 ( 0.06 7.00 7.59

expt

60 MPa

44

MD

5.78 ( 0.07 6.45

7.88 8.52 6.41

expt

44

MD

5.15 ( 0.07 5.73

6.73 6.82

5.16 7.01

TABLE 4: Isothermal Compressibility, KT, of [C4mim+][Tf2N-] and [C6mim+][Tf2N-] at Atmospheric Pressure from Experimental Measurements and NPT MD Simulations κT (×104 MPa-1)

ionic liquid 298.15 K

308.15 K

318.15 K

44

44

44

expt [C4mim+][Tf2N-] [C6mim+][Tf2N-]

5.26 ( 0.08 5.61 ( 0.03

MD 3.35 4.16

expt

5.54 ( 0.02 5.57 ( 0.01

from DFT calculations for the [C6mim+] cation and the [Tf2N-] anion is shown in Figure 2. Results for [C4mim+] and [C8mim+] cations are very similar to those for [C6mim+] and are not shown here. Blue areas in the cation represent the most positive electrostatic regions, while the red color in the anion displays the most negative ones. Ring atoms in between the two nitrogen atoms (C2 and H1) are more electropositive than the other carbon and hydrogen atoms in the ring (C4, C5, H2, and H3). Atoms on the alkyl chain are considerably less positive that those in the ring. On the other hand, the most negative areas in the anion are concentrated in the oxygen and nitrogen atoms. Hence, it is expected that the negative regions of the anion are placed close to the positive areas of the cation. Keeping these observations in mind, various different conformers were found for the different [Cnmim+][Tf2N-] ion pairs by matching the most negative and positive regions for the anion-cation pairs. 3.1.2. Geometries and energies for isolated [Cnmim+][Tf2N-] ion pairs. In all conformers examined, the energy of formation for an ion pair, defined as

∆Eion pair ) E[Cnmim+][Tf2N-] - (E[Cnmim+] + E[Tf2N-])

(2)

is highly exothermic, with ∆Eion pair > 300 kJ/mol. Relative energies for the most stable conformers of the ion pairs are shown in Figure 3. The values in parentheses correspond to calculations where the zero-point energy corrections were taken into account. The differences between the two sets of calculations are small. Conformer III was set arbitrarily as the reference conformer. In addition, the most important geometrical parameters of these conformers are presented in Table 1. The conformers are characterized by two different geometrical parameters, the C2-N4 distance and the dihedral angle X-H1-N4-S2, where X is the centroid of the imidazolium ring. A value of this angle around 180° indicates a perpendicular orientation of the anion with respect to the imidazolium ring, whereas values close to 0° denote a parallel orientation. In structures I-III and V, the nitrogen atom N4 of the anion is pointing out to the C2 carbon of the imidazolium ring (see C2-N4 distances in Table 1). The difference between these conformers is that in structures I-III, the anion is placed in an almost perpendicular plane with respect to the imidazolium ring, whereas in the V conformer, the anion is in a roughly parallel orientation with respect to the imidazolium ring (see the X-H1-N4-S2 dihedral angle in Table 1). On the contrary, conformer IV has the O2 and O4 oxygen atoms pointing out to

MD

expt

4.78 4.17

5.83 ( 0.03 5.65 ( 0.02

328.15/333.15 K MD 4.70 4.54

expt44

MD

5.87 ( 0.04

4.93 4.84

the imidazolium ring (see C2-N4, C2-O2, and C2-O4 distances in Table 1). This conformer is the lowest-energy one by around 3.1-7.2 kJ/mol with respect to the other conformers. On the basis of the DFT calculations, it can be concluded that the alkyl chain length of the cation does not affect the relative stability between the different conformers. Structures for [C4mim+][BF4-] and [C4mim+][I-] ILs in which the anion is positioned on top of the cation have been experimentally proposed.41 We have found that the VI conformers of the [Cnmim+][Tf2N-] ion pairs place the anion on top of the cation with a higher energy with respect to conformer IV by around 3.0 kJ/mol. In accordance with these energetic differences, it is expected that the conformer IV would prevail with respect to the others in the bulk liquid. Consequently, the lowest-energy conformer IV has been used in the construction of the initial configurations of the ILs in the bulk phase that were used for the MD simulations. The averaged electrostatic charges over all ion pair conformers calculated by the DFT/CHelpG method were used in the development of the Coulombic potential parameters in the force field. The electrostatic charges for all conformers are available in Tables S2-S4 in the Supporting Information. Other researchers have used the atomic charges by taking into account the isolated geometries for the anion and the cation.42 Experimentally, two peaks have been observed in the X-ray photoelectron spectrum (XPS) of [C4mim+][Tf2N-], which were assigned to the two nitrogen atoms (N1 and N3 of Figure 1) in the imidazolium ring and to the nitrogen in the sulfoimide anion (N4 in Figure 1).43 Caporali et al.43 attributed these two different peaks to the presence of a negative charge localized on the N4 atom due to the sulfur atoms. The electrostatic charges calculated in this work for the N1, N3, and N4 atoms are 0.153, 0.063, and -0.781 electrons (averaged over all conformers), which are in a very good agreement with the experimental results. Furthermore, XPS experiments showed just one set of sulfur atoms, which is in close agreement with the similar electrostatic charges calculated for S1 and S2 here of about 1.166 and 1.124 electrons, respectively. For the carbon atoms, it is clear both experimentally and theoretically that there are two sets of atoms, the first one of very electropositive carbons linked to fluoride atoms (C41 and C42) and the second one of less electropositive or even electronegative carbon atoms belonging to the cation (see electrostatic charges in the tables in the Supporting Information). Similar conclusions can be drawn from the electrostatic charge

7216


Logotheti et al.

Figure 6. Center of mass (COM) radial distribution functions for the three ILs studied at 318.15 K and 0.1 MPa; (a) anion-anion, (b) cation-anion, and (c) cation-cation.

Figure 5. Probability distribution of various dihedral angles in (a and b) the cation and (c) the anion of [C6mim+] [TF2N-] at 298.15 K and 0.1 MPa from MD simulations.

analyses for the [C6mim+][Tf2N-] and [C8mim+][Tf2N-] ILs for which no experimental results are available. This good

agreement between experimental data and calculations is a clear validation of the approach used here for the parametrization of the force field in terms of atom charges. 3.2. Molecular Dynamics Simulations. 3.2.1. Volumetric and Thermodynamic Properties of the ILs. In Table 2, mass density predictions from the NPT MD simulations for the three ILs are shown. For a given temperature and pressure value, an increase in the alkyl tail size of the cation results in a decrease in density. In Figure 4, MD predictions are compared against experimental data from the literature.44,45 The agreement is very good in all cases; MD predictions are within less than 1% of the experimental data, while in most cases, the deviation is on the order of 0.2-0.5%. It should be emphasized here that this very good agreement between experiments and simulation was achieved without any force field parameter adjustment to the experimental data. For comparison, similar density calculations using the force field of Carnogia Lopes and Pa´dua19 result in deviations from experimental data on the order of 3%.19,21 An agreement with experimental data similar to the one reported here was only obtained by an arbitrary refit of the Lennard-Jones parameters of the model.21



Figure 7. Atom-atom RDF for [Cnmim+][Tf2N-] ion pairs (n ) 4, 6, and 8) at 318.15 K and 0.1 MPa. Molecular structures correspond to the isolated conformers optimized at the quantum DFT level.

Figure 8. C2-N4 radial distribution function for [Cnmim+][Tf2N-] ion pairs (n ) 4, 6, and 8) at 318.15 K and 0.1 MPa. The left side structure represents approximately the first coordination sphere (molecules in ball-stick style) of a randomly chosen ion pair (VdW model). The colors of the arrows and ball-stick models match the distances in the figure, indicating intermolecular interactions. The right side figure highlights the intermolecular interactions of different ion pairs raising a like-ordered structure of anions (red) and cations (blue) in the bulk liquid.

(4)

3.2.2. Conformational Properties of ILs. A major issue regarding the conformational properties of IL molecules arises from the fact that the electrostatic potential has been parametrized based on DFT calculations while all other force field parameters were taken from the literature. In this respect, a thorough analysis of the distribution of some major dihedral angles is necessary. In Figure 5, the probability distributions of N3-C7-C8-C9, C7-C8-C9-C10, and C2-N3-C7-C8 dihedral angles for [C6mim+] and of the S1-N4-S2-C42 dihedral angle for [Tf2N-] at 298.15 K are shown.

In eqs 3 and 4, 〈 〉 denotes statistical ensemble averages calculated during simulation. In Tables 3 and 4, MD predictions and experimental data44,45 are shown for aP and κT of the three ILs at different pressures and temperatures. The agreement between simulation and experiments is reasonable in all cases, verifying that the force field developed based on quantum DFT calculations can be used to predict even the derivative bulk properties of the ILs.

The distribution of the N3-C7-C8-C9 angle exhibits a maximum for the trans (or anti) conformation at 180° and two local maxima for the two gauche conformers at -60 and 60°. This is in agreement with the DFT calculations at the B3LYP/ 6-311+G* level that predicted three minima for the torsion energy profile at 180, -60, and 60°, respectively (not shown here). The trans conformer is more stable by about 2-3 kJ/mol than the gauche ones, thus resulting in a more populated trans

MD data were also used to calculate additional thermodynamic properties such as the isobaric thermal expansion coefficient, aP, through the expression

aP )

( )

1 ∂V 1 ∆〈V〉 ) V ∂T P 〈V〉 ∆T

( )

(3)

P

and the isothermal compressibility, κT, through the expression

κT ) -

( )

1 ∂V 1 ∆〈V〉 )V ∂P T 〈V〉 ∆P

( )

T

7218


Logotheti et al.

Figure 9. RDFs between two carbon atoms (C2 and C4) in the imidazolium ring of [C4mim+] and various atoms (N4, S1, and O2) in the [Tf2N-] anion at 318.15 K and 0.1 MPa.

Figure 10. Radial distribution functions for the carbon atoms of the ring of the [C4mim+] and an oxygen atom of the [Tf2N-] at 318.15 K and 0.1 MPa.

conformation. The fraction of the trans conformation as calculated from Figure 5 is xA ) 0.78 (using the nomenclature from ref 20), slightly larger than the experimental value of xA ) 0.65 based on Raman experiments.46 Canongia Lopes and Pa´dua, on the basis of MD simulations,20 and Berg et al.,46 on the basis of quantum mechanics calculations, predicted a lower value for xA, equal to 0.39 and 0.50, respectively. The difference between the two sets of MD data can be attributed, among others, to the geometry optimization method used in each case. The probability distribution of the C7-C8-C9-C10 dihedral angle (Figure 5a) is similar to that of the N3-C7-C8-C9 angle, except for an increase in the trans population for the former angle. The probability distribution of the C2-N3-C7-C8 angle shown in Figure 5b is in good agreement with those reported by Canongia Lopes and Pa´dua.20 Finally, the S1-N4-S2-C42 dihedral angle distribution exhibits two broad peaks centered at 95 and 230° (Figure 5c). This distribution agrees with the MP2/cc-pVTZ(-f) calculations previously reported by Canongia Lopes and Pa´dua.19 In addition, the peak at 95° is slightly higher than the second one, which is in agreement with an energy difference of 4-5 kJ/mol between the two energy minima as reported by the same authors. 3.2.3. Liquid Bulk Structure. A detailed analysis of the microscopic bulk structure of the ILs is based on the radial distribution function (RDF). In Figure 6, the center of mass RDFs (COM-RDF) are shown for the (a) anion-anion, (b) cation-anion, and (c) cation-cation of the three ILs at 318.15 K and 0.1 MPa. In all cases, the cation-anion COMRDF exhibits a strong first peak at 4.9-5.3 Å, indicating a preferential cation-anion interaction over anion-anion and cation-cation interactions. The cation-anion COM-RDF

Figure 11. Second-order autocorrelation function P2(t) for various bond vectors of (top) the cation in [C8mim+][Tf2N-] and (bottom) the anion in [C4mim+][Tf2N-] at 318.15 K and 0.1 MPa.

exhibits also a second less intense but broader peak at around 12 Å, while ordering persists at longer distances close to or even beyond 20 Å, which is due to the long-range Coulombic interactions. This behavior has been observed for other ILs as well, both in experiments47-50 and in simulations.9,39,40,51-54 Anion-anion COM-RDFs exhibit a split first peak at 7.8 and 10.3 Å. A second peak for this pair is observed at 16.4 Å. Finally, the cation-cation COM-RDFs exhibit a broad first peak at 8.5-9.5 Å with a low intensity compared to that of the other two pair types. In all COM-RDFs shown in Figure 6, the intensity of a given peak decreases, and its width increases as the alkyl tail size of the cation increases (from butyl to hexyl to octyl). This is likely due to the higher conformational flexibility of the longer alkyl chains resulting in a delocalization of the COM of the cation. At the macroscopic level, this behavior is reflected in lower mass densities for ILs with longer alkyl tails (see Table 2). The alkyl tail size effect on the calculations in Figure 6 is in good agreement with MD simulation results for [Cnmim+][PF6-] ILs where n ) 4, 6, 8, and 12.23 RDF data can be used to calculate the coordination number N for a given ion through the expression


N ) 4πF

∫0rshell r2g(r)dr


(5)

where F is the bulk number density and rshell is the position of the first minimum of the RDF.9,11,13 The results are very similar for all ILs studied; N turns out to be 1, and the ion is of opposite sign compared to the central one at about 5.3 Å. Using rshell equal to 9 Å (which corresponds to the first minimum of the RDF), the calculated coordination number for the cation-anion first solvation shell ranges from 5.4 to 5.65. The smaller value of 5.45 corresponds to [C8mim+][Tf2N-], which agrees with the larger size of the cation and the more diffusive behavior that exhibits. The same conclusion was reached by other researchers for pyridinium-based ILs.39 Further insight into the bulk liquid structure can be extracted from the atom-atom RDFs. RDFs for the nitrogen atom in the anions and two carbon atoms in the imidazolium ring of the cations are shown in Figure 7 at 318.15 K and 0.1 MPa. Four well-defined peaks are identified in the lefthand side figure (C2-N4 RDF) at 3.3, 4.8, 7.3, and 11.7 Å. For the right-hand side figure (C4-N4 RDF), an intense peak is observed at 5.4 Å associated with a shoulder at 6.7 Å. In addition, a broad peak is found at 12.5-13.5 Å. The first peaks in each atom-atom RDF (at 3.3 Å in the left-hand

Figure 12. mKWW fit to P2(t) simulation data for the C7-C8 bond of the [C6mim+][Tf2N-] IL at 318.15 K and 0.1 MPa.

Figure 13. Time dependence of the P2(t) orientation autocorrelation function for various bond vectors of [C4mim+][Tf2N-] at four different temperatures. Black lines indicate data at 298.15 K, red lines data at 308.15 K, green lines data at 318.15 K, and blue lines data at 328.15 K.

Figure 14. Torsion autocorrelation function Pφ(t) for various torsion angles in (top) the cation of [C8mim+][Tf2N-] and (bottom) the anion of [C4mim+][Tf2N-] at 318.15 K and 0.1 MPa.

side and 5.4 Å in the right-hand side figure) can be assigned to intramolecular C2-N4 and C4-N4 interactions, which are very similar to those found by DFT calculations. Furthermore, the peaks corresponding to lowest-energy conformer IV are the most intense. Hence, a good match between the isolated ion pair structure calculated by DFT and the structure of the ion pair in the bulk liquid is obtained. Nevertheless, there are two peaks at 7.3 and 11.7 Å on the left-hand side and one peak at 12.5-13.5 Å on the right-hand side of Figure 7 that cannot be assigned to intramolecular interactions only and should be primarily attributed to intermolecular interactions. The effect of the alkyl tail on the intensity of the primary peaks in Figure 7 is in agreement with MD simulations for [Cnmim+][PF6-] ILs where n ) 6, 8, 10, and 12.26 In Figure 8, the first-neighbor intermolecular interactions between the different ion pairs in the bulk liquid are shown. A close examination of this figure allows one to assign the two peaks at 7.3 and 11.7 Å (left picture in Figure 8). Furthermore, neighboring cation-anion pairs are structured in layers, as shown in the right picture of Figure 8. Here, the anions are represented in red and the cations in blue; therfore, one can clearly identify the anion-cation intermolecular interactions. This structure is manifested in the shoulder of the first peak at 5.6 Å. However, for the [C8mim+][Tf2N-] IL, this shoulder is not observed, likely due to the increase of the alkyl tail length that makes this structure less likely. Clearly, additional work is needed in order to clarify further these first observations. In Figure 9, the RDFs between each of the two carbons C2 and C4 of the imidazolium ring of [C4mim+] and the N4, S1, and O2 atoms of [Tf2N-] at 318.15 K and 0.1 MPa are shown. In both cases, the O2 atom is closer to the carbon atoms, followed by the S1 and then by the N4 atom. The O2 is covalently bonded to a single atom only (S1), and at the same time, it possesses the largest negative charge on the anion; thus, it gets closer to the carbon atoms of the ring that possess the largest positive charge in the cation. At the same time, the S1 atom also gets close to the C2 and C4 atoms. The C2-N4 and C4-N4

7220


Logotheti et al. TABLE 5: Mean Decorrelation Time, τC, for Various Bonds of [C4mim+][Tf2N-], [C6mim+][Tf2N-], and [C8mim+][Tf2N-] Calculated from the mKWW Fit to the P2(t) Dataa bond

τC (ns) 298.15 K

Figure 15. Time dependence of the Pφ(t) function for the torsion angle C41-S1-N4-S2 (squares) and F1-C41-S1-O2 (circles) of the [C4mim+][Tf2N-] IL at 318.15 K. The red curve corresponds to the modified stretched exponential (mKWW) and the green curve to the stretched exponential (KWW) functions fitted to simulation results.

N1-N3 N1-C6 N3-C7 C7-C8 C9-C10 C41-F2 C41-S1 S1-O2 S1-N4

[C4mim >200b 4.4 ( 0.3 >300b 2.5 ( 0.3 >50b 5.9 ( 0.5 1.8 ( 0.9 0.61 ( 0.13 1.1 ( 0.4 0.38 ( 0.06 0.74 ( 0.09 0.409 ( 0.005 11.8 ( 8.3 1.11 ( 0.08 11.1 ( 4.2 1.35 ( 0.16 16.4 ( 13.3 1.32 ( 0.13 298.15 K

N1-N3 400b N1-C6 >1000b N3-C7 20.7 ( 4.5 C7-C8 4.4 ( 2.1 C9-C10 4.7 ( 2.5 C11-C12 1.5 ( 0.6 C41-F2 0.76 ( 0.11 C41-S1 3.7 ( 1.6 S1-O2 3.9 ( 1.3 S1-N4 3.6 ( 1.4 298.15 K N1-N3 N1-C6 N3-C7 C7-C8 C9-C10 C11-C12 C13-C14 C41-F3 C41-S1 S1-O2 S1-N4

Figure 16. Torsion autocorrelation function Pφ(t) for various torsion angles in the cation (top) and anion (bottom) of the [C6mim+][Tf2N-] IL at 298.15 (black), 308.15 (red), 318.15 (green), and 333.15 K (blue).

RDFs exhibit a first peak at relatively longer distances due to the fact that N4 is shielded by the two Tf groups. Similar behavior has also been observed by other researchers in the simulation of pyridinium-based ILs.39 In Figure 10, the RDFs of the three carbon atoms in the imidazolium ring (C2, C4, and C5) with O2 in the anion at 318.15 K are shown. Clearly, O2 is located preferentially closer to C2 than in any of the other two carbon atoms. This is in full agreement with the characteristic low-energy conformers identified by the DFT calculations. RDFs can be also used to elucidate the heterogeneity of ILs and primarily the formation of nanostructures. Preliminary calculations of the RDF between terminal carbon atoms in the alkyl tails of the cations based on the MD simulations confirmed such a phenomenon. Further detailed analysis of the effect of temperature, pressure, and alkyl tail size on this phenomenon will be presented in a future work. 3.2.4. Local Dynamics of IL Molecules. The local (segmental) dynamics of the IL molecules was examined through the investigation of bond decorrelation and torsion angle decorrelation. The

308.15 K

17.6 ( 6.2 9.8 ( 3.1 20.0 ( 8.7 2.97 ( 0.09 9.9 ( 4.7 10.0 ( 5.3 3.4 ( 1.1 0.72 ( 0.08 3.73 ( 0.97 6.4 ( 3.4 4.4 ( 1.2

318.15 K +

328.15 K

-

][Tf2N ] 2.3 ( 0.3 0.92 ( 0.49 1.34 ( 0.08 0.87 ( 0.11 2.9 ( 0.7 2.5 ( 0.6 0.31 ( 0.04 0.21 ( 0.03 0.20 ( 0.01 0.13 ( 0.02 0.26 ( 0.02 0.186 ( 0.002 0.63 ( 0.06 0.39 ( 0.04 0.68 ( 0.13 0.43 ( 0.04 0.72 ( 0.09 0.48 ( 0.03

308.15 K 318.15 K 333.15 K [C6mim+][Tf2N-] >200b 2.6 ( 0.2 1.54 ( 0.01 >400b 1.42 ( 0.09 0.829 ( 0.006 13.0 ( 1.7 3.3 ( 0.4 1.73 ( 0.06 1.3 ( 0.2 0.70 ( 0.08 0.38 ( 0.03 1.5 ( 0.1 0.67 ( 0.08 0.37 ( 0.04 0.59 ( 0.03 0.26 ( 0.03 0.15 ( 0.01 0.44 ( 0.02 0.27 ( 0.01 0.416 ( 0.006 1.4 ( 0.3 0.63 ( 0.02 0.61 ( 0.05 1.5 ( 0.3 0.74 ( 0.06 0.43 ( 0.06 1.6 ( 0.3 0.76 ( 0.02 0.46 ( 0.04 308.15 K 318.15 K 328.15 K [C8mim+][Tf2N-] 15.0 ( 8.2 5.5 ( 1.3 2.7 ( 0.6 8.0 ( 3.1 3.0 ( 0.6 1.6 ( 0.2 13.7 ( 6.4 6.2 ( 2.2 2.9 ( 0.6 1.9 ( 0.4 1.7 ( 0.3 1.02 ( 0.04 5.3 ( 1.5 3.0 ( 1.1 1.47 ( 0.11 3.9 ( 1.4 2.3 ( 0.5 1.16 ( 0.07 1.7 ( 0.5 1.0 ( 0.3 0.56 ( 0.07 0.57 ( 0.08 0.34 ( 0.04 0.238 ( 0.008 3.0 ( 1.4 1.2 ( 0.4 0.63 ( 0.08 3.8 ( 1.7 1.5 ( 0.3 0.73 ( 0.06 2.6 ( 0.9 1.3 ( 0.3 0.75 ( 0.11

a The statistical uncertainty is based on block average analysis. All data refer to 0.1 MPa. b The value is too large to be estimated accurately.

reorientation of bonds in the cations and anions can be quantified through a second-order autocorrelation function of the form

1 P2(t) ) {3〈[uCH(t) · uCH(0)]2 〉 - 1} 2

(6)

where uCH(t) is the unit vector along a bond at time t. In Figure 11, P2(t) as a function of time is shown for various bonds of a cation ([C8mim+]) and an anion at 318.15 K and 0.1 MPa. Similar behavior was observed for the other cations as well. Two different types of bonds in the cation can be identified based on how fast P2(t) approaches 0. Bonds between atoms in the imidazolium ring (N1-N3, N1-C6, and N3-C7) are relatively stiff and decorrelate slow. On the other hand, bonds between carbon atoms in the alkyl tail (C7-C8, C9-C10, C11-C12, and C13-C14) are significantly more mobile and decorrelate fast. This is partly due to the fact that the partial charges on these carbon atoms are smaller than the charges in the ring atoms. The fastest decorrelating bond is the one connecting the last two carbon atoms in the tail. In the anion, the terminal bond C41-F2 decorrelates faster because it can move freely and the partial charges are relatively small. At the same time, the S1-N4 bond is stiffer since it is an internal bond in the anion and partial charges in the two atoms are relatively higher.



TABLE 6: Mean Decorrelation Time, τC, for Various Torsion Angles of [C4mim+][Tf2N-], [C6mim+][Tf2N-], and [C8mim+][Tf2N-] Calculated from the mKWW or KWW Fit to the PO(t) Dataa torsion angle

τC (ns) 298.15 K

308.15 K

318.15 K +

C2-N3-C7-C8 N3-C7-C8-C9 C7-C8-C9-C10 C41-S1-N4- S2 F1-C41-S1-O2 F1-C41-S1-N4

[C4mim ][Tf2N ] 0.11 ( 0.02 0.064 ( 0.005 0.17 ( 0.03 0.097 ( 0.013 0.068 ( 0.009 0.059 ( 0.012 0.26 ( 0.06 0.18 ( 0.01 1.27 ( 0.07 0.82 ( 0.01 1.3 ( 0.2 0.92 ( 0.10

0.38 ( 0.16 0.25 ( 0.03 0.11 ( 0.01 1.2 ( 0.4 2.1 ( 0.2 2.1 ( 0.2 298.15 K

C2-N3-C7-C8 N3-C7-C8-C9 C7-C8-C9-C10 C8-C9-C10-C11 C9-C10-C11-C12 C41-S1-N4-S2 F3-C41-S1-O2 F3-C41-S1-N4

308.15 K

318.15 K [C6mim+][Tf2N-] 0.24 ( 0.02 0.11 ( 0.02 0.31 ( 0.08 0.19 ( 0.01 0.15 ( 0.02 0.088 ( 0.012 0.11 ( 0.02 0.078 ( 0.009 0.068 ( 0.015 0.044 ( 0.003 0.24 ( 0.06 0.13 ( 0.05 1.4 ( 0.3 0.88 ( 0.05 1.2 ( 0.2 0.84 ( 0.07

0.51 (0.14 0.74 ( 0.52 0.21 ( 0.08 0.18 ( 0.03 0.090 ( 0.012 0.67 ( 0.25 1.48 ( 0.08 1.9 ( 0.3 298.15 K

C2-N3-C7-C8 N3-C7-C8-C9 C7-C8-C9-C10 C8-C9-C10-C11 C9-C10-C11-C12 C10-C11-C12-C13 C11-C12-C13-C14 C41-S1-N4-S2 F2-C41-S1-O2 F2-C41-S41-N4 a

328.15 K

-

308.15 K

318.15 K [C8mim+][Tf2N-] 0.34 ( 0.03 0.26 ( 0.01 0.78 ( 0.32 0.36 ( 0.04 0.21 ( 0.01 0.18 ( 0.05 0.49 ( 0.10 0.25 ( 0.04 0.36 ( 0.18 0.13 ( 0.03 0.080 ( 0.022 0.078 ( 0.006 0.074 ( 0.040 0.065 ( 0.020 1.23 ( 0.80 0.53 ( 0.23 1.28 ( 0.07 0.87 ( 0.12 1.32 ( 0.05 0.87 ( 0.08

0.67 ( 0.26 1.7 ( 0.8 0.39 ( 0.13 0.65 ( 0.15 0.42 ( 0.26 0.17 ( 0.05 0.067 ( 0.011 1.22 ( 0.59 1.33 ( 0.17 1.74 ( 0.35

0.048 ( 0.007 0.065 ( 0.018 0.040 ( 0.004 0.07 ( 0.01 0.681 ( 0.03 0.67 ( 0.03 333.15 K 0.056 ( 0.007 0.131 ( 0.001 0.054 ( 0.003 0.054 ( 0.005 0.032 ( 0.001 0.094 ( 0.007 0.57 ( 0.03 0.55 ( 0.01 328.15 K 0.185 ( 0.004 0.29 ( 0.05 0.12 ( 0.01 0.18 ( 0.01 0.087 ( 0.005 0.049 ( 0.004 0.032 ( 0.002 0.28 ( 0.09 0.73 ( 0.07 0.66 ( 0.01

The statistical uncertainty is based on block average analysis. All data refer to 0.1 MPa.

The data for P2(t) shown in Figure 11 can be described accurately with the modified Kohlrausch-Williams-Watts (mKWW) function of the form55

[ ( )]

t PmKWW(t) ) R exp τ0

[ ( )]

+ (1 - R) exp -

t

τC )

β

τKWW

(7) This mKWW function consists of two terms. The first one describes a fast exponential decay with amplitude R associated with small perturbation of torsion angle values around skeletal bonds and with the bond and bond-angle bending vibrations of skeletal and pendent bonds near their equilibrium values, with characteristic time τ0. The second term is a slower stretched exponential (KWW) decay associated with cooperative conformational transitions in the ion, with τKWW being the characteristic correlation time and β the stretching exponent. In Figure 12, a typical fit of mKWW function to raw P2(t) simulation data is shown for the C7-C8 bond of [C6mim+][Tf2N-] at 318.15 K and 0.1 MPa. The agreement is excellent over the entire time range examined. A temperature increase results in faster bond decorrelation. In Figure 13, P2(t) for N3-C7 for the cation and C41-F2 for the anion of [C4mim+][Tf2N-] at four different temperatures is shown. Similar behavior was observed for the other bonds and all ILs studied. The area under P2(t) defines the correlation time τC for the reorientational motion of the bond studied

τC )

∫0∞ P2(t)dt

relaxation of the corresponding bond can be derived analytically using eqs 7 and 8, so that

(8)

Given that the mKWW function, PmKWW(t), provides an excellent fit to P2(t) data, the mean decorrelation time of the orientational

∫0∞ P2(t)dt)∫0∞ PmKWW(t)dt) 1 1 Rτ0 + (1 - R)τKWW Γ β β

()

(9)

In this way, simulation data were used to evaluate τC for the various bonds at different conditions, and results are reported in Table 5. Clearly, the decorrelation time varies by almost 3 orders of magnitude between bonds in the imidazolium ring (primarily N1-N3 and N1-C6) of the [C4mim+] and [C6mim+] cations and bonds in the end of the alkyl tail. For the very slow decorrelating bonds, much longer runs are needed in order to obtain accurate estimates of τC. The significant variation in the decorrelation of bonds depending on their position in the IL molecule resembles the behavior in polymers where backbone bonds are significantly stiffer than terminal bonds in short branches, as evaluated from both simulations55,56 and experiments.55,57 Interestingly enough, the bonds in the imidazolium ring of [C8mim+] are significantly less stiff,with a characteristic decorrelation time on the same order of magnitude as the decorrelation time for bonds in the alkyl tail. Experimental data27 indicated that the increase of alkyl tail size in cations has a nonmonotonic effect on the melting temperature (Tm) of the IL. For imidazolium-based ILs, Tm is on the order of 300 K, whereas the exact value depends on the size of the alkyl tail and on the anion.27 One can argue here that the decrease of decorrelation time for [C8mim+] at 298 K compared to that for [C4mim+] and [C6mim+] is an indication that this IL has a lower Tm compared to the other two.

7222


Logotheti et al.

Figure 18. Self-diffusion coefficient D of [C4mim+][Tf2N-], [C6mim+][Tf2N-], and [C8mim+][Tf2N-]. (Top) Cation diffusion coefficients; (bottom) anion diffusion coefficients. In all cases, the lines are experimental data,16 and the points are MD calculations.

typically to model data for long chain molecules where synergistic motion is important.55 On the other hand, Pφ(t) data for torsion angles in the anions exhibit different behavior, as can be seen in Figures 14 and 15. The torsion angle C41-S1-N4-S2 that consists of internal atoms to the anion exhibits also a mKWW dependence with time, whereas torsion angles F1-C41-S1-O2 and F1-C41-S1-N4 that contain terminal atoms can be described with a simple KohlrauschWilliams-Watts (KWW) function of the form

[ ( )]

PKWW(t) ) (1 - R) exp -

Figure 17. Running average of β defined in eq 13 for the case of [C4mim+][Tf2N-] at 298.15 and 308.15 K and of [C8mim+][Tf2N-] at 298.15 K.

An additional approach to quantify the segmental dynamics of IL molecules is based on the torsional dynamics that is quantified through the torsion autocorrelation function, defined as

Pφ(t) )

〈cos(φ(t))cos(φ(0))〉 - 〈cos(φ(0))〉2 〈cos(φ(0))cos(φ(0))〉 - 〈cos(φ(0))〉2

(10)

where φ(t) is a torsion angle at time t. In Figure 14, Pφ(t) simulation data are shown for the torsion angles of the [C8mim+] cation and the [Tf2N-] anion examined at 318.15 K and 0.1 MPa. Torsion angles with respect to atoms in the alkyl tail decorrelate faster than torsion angles of atoms near the imidazolium ring, in agreement with the previous analysis for the bond decorrelation. The torsion autocorrelation function data for the cations can be described accurately by the mKWW function, which is used

t τKWW

β

(11)

Representative fits of these functions to simulation data are shown in Figure 15. The Pφ(t) data for the anions in all three ILs studied display the same behavior. The effect of temperature on Pφ(t)is shown in Figure 16 for [C8mim+][Tf2N-]. As temperature increases, torsion angles decorrelate faster. This behavior is similar for the other two ILs as well. The time integral of Pφ(t) (eq 8) provides the correlation time τC for the torsional reorientational motion. Given that the mKWW function PmKWW(t) provides an excellent fit to Pφ(t) data for all cation and some anion torsion angles, the mean correlation times for these torsion angles can be obtained from eq 9. For those torsion angles that follow a KWW behavior, an analogous expression for the calculation of τC is used. The calculated correlation times and their statistical uncertainty at all temperatures for the ILs examined are shown in Table 6. Overall, a comparison of the decorrelation times reported in Tables 5 and 6 indicate that torsion angles, on average, decorrelate faster than bonds. 3.2.5. Self-Diffusion Coefficient of ILs. Long MD runs allow reliable calculation of the self-diffusion coefficient of the cation and of the anion of an IL based on the mean square displacement (MSD) of each ion in the Fickian regime using Einstein’s equation


1d ri(0)|2 〉 〈|b r (t) - b 6 dt i


4. Conclusions

in all cases. The microscopic structure was studied thoroughly, primarily in terms of radial distribution functions and in terms of the coordination number. It was shown that the change of the alkyl tail size, from butyl to hexyl to octyl, has relatively little effect on the radial distribution functions. In addition, the microscopic structure in bulk liquid has a lot of similarities with minimum-energy conformers identified for isolated ion pairs. This provides a nice validation of the approach used to incorporate the density functional theory results into the molecular model. Of course, a major concern here refers to the omission of polarization and many body interactions. Explicit inclusion of such interactions in the calculations would result in a significant increase in the computational time and would make the DFT and long MD simulations performed here intractable. Although it is difficult to quantify accurately the effect of these higher-order interactions into the bulk properties examined here, one may argue that the effective potential developed in this work provides an accurate prediction of structure, thermodynamic, and dynamic properties of ILs against experimental data, and thus, the contribution of higher-order terms is expected to be relatively small. Thermodynamic properties that include mass density, the isobaric thermal expansion coefficient, and isothermal compressibility were calculated and found to be in good to very good agreement with literature experimental data. In the case of mass density, the increase of the alkyl tail size results in reduction of the mass density for a given temperature and pressure. The new force field is an improvement over an earlier one in terms of accuracy in the prediction of IL thermodynamic properties.19 The local segmental dynamics of IL molecules was evaluated based on bond decorrelation and torsion angle decorrelation. In terms of the former calculation, it was found that bond decorrelation for the case of [C4mim+] and [C6mim+] varies by almost 3 orders of magnitude from the very stiff bonds associated with the imidazolium ring to the very mobile bonds at the end of the alkyl tail. Such behavior resembles the variability in bond decorrelation in long branched polymer chains. On the other hand, decorrelation of various types of torsion angles was found to exhibit a significantly smaller distribution in terms of the characteristic decorrelation time, τC. Finally, long MD calculations provided accurate prediction of the cation and anion self-diffusion coefficients over the entire temperature range examined. For both local and global dynamics, it became obvious based on relaxation times and selfdiffusion coefficients that as the alkyl chain of the cation increases, the dynamics slows down due to the tail aggregation effect. Future work will focus on spatial heterogeneities in bulk ILs and the effect of pressure on the local segmental dynamics and the self-diffusion coefficients.

In this work, the molecular structure, the thermodynamic properties, and the segmental dynamics of three ILs consisted of the same anion ([Tf2N-]) and different cations that contain an imidazolium ring and alkyl tails of variable size ([C4mim+], [C6mim+] and [C8mim+]) were investigated using various molecular modeling techniques. Quantum calculations using density functional theory were performed on isolated ion pairs, and a series of low-energy conformers were identified. In parallel, the electron density distribution on the ions was accurately calculated and used in the molecular force field for the specific ILs. Subsequently, long MD simulations were performed for bulk ILs at various temperatures and pressures. The CHARMM force field for the bonded and Lennard-Jones interactions was used

Acknowledgment. Financial support provided by INTAS through Research Project No. 05-1000008-8020 on Development of Sustainable Industrial Processes: Experimental, Theoretical and Computational Investigation of Thermodynamic Properties and Phase Equilibria of Ionic Liquid Mixtures is gratefully acknowledged. J.R. thanks CSIC and CICYT for financial support through an I3P tenure track position and MAT20060400 project, respectively. The authors also acknowledge Centro Tećnico de Informa´tica (CTI-CSIC) and Centro de Investigaciones Energe´ticas, Medioambientales y Tecnolo´gicas (CIEMAT), both in Madrid, Spain, for the use of their computational resources for DFT and MD calculations. Finally, we are grateful to Dr. Loukas Peristeras of Scienomics SARL in Paris, France, for the technical assistance provided regarding the

D ) lim

tf∞

(12)

In order to ensure that the system’s diffusion is in the Fickian regime, Maginn and co-workers39 proposed that the following quantity is calculated from the MD data

β(t) )

d log〈|b r i(t) - b r i(0)|2 〉 d log(t)

(13)

where β ) 1 in the Fickian regime. At short times, ballistic motion dominates, and therefore, β ) 2, while at intermediate time scales, ILs exhibit glass-like behavior, and therefore, β < 1. On the basis of relatively short MD simulations (on the order of 5 ns) for pyridinium-based [Tf2N-] ILs, Maginn and co-workers concluded that the ILs exhibited glassy dynamics at least at low temperatures, and only at 423 K β approached 1 after approximately 2 ns. In Figure 17, representative results concerning the running average of β for [C4mim+][Tf2N-] at 298.15 and 308.15 K and for [C8mim+][Tf2N-] at 298.15 K are shown for the first 10 ns of MD simulation. Clearly, the normal diffusive regime is attained within a few nanoseconds, ensuring the reliability of diffusion coefficient values obtained from simulation runs on the order of 50 ns. In Figure 18, experimental data16 and MD predictions are shown for the self-diffusion coefficient of the various cations and anions examined in this work as a function of temperature. The increase of alkyl tail size results in the decrease of the diffusion coefficient of the cation and also that of the corresponding anion, but to a slightly less extent. The reduction in the diffusion coefficient of cations with the increase of alkyl tail size has been linked to the spatial heterogeneity of ILs, which increases for larger alkyl tails.58 In all cases, MD calculations are in very good agreement with experimental data. Some deviation at 298.15 K for the diffusion coefficient of [Tf2N-] can be improved with longer MD runs. The percentage average absolute deviation between experimental data and MD calculations over all D values is 13%, of the same magnitude as the estimated statistical uncertainty of the calculations. The accuracy of the force field proposed here in predicting D compares very favorably with previous models. Ko¨ddermann et al.21 showed that D values for [C2mim+][Tf2N-] in the temperature range of 273-373 K calculated using the force field of Canongia Lopes et al.14 differ by up to an order of magnitude from experimental data. A refinement of the model based on the density and diffusion coefficient of [C2mim+][Tf2N-] resulted in D values for [C1mim+][Tf2N-] up to [C8mim+][Tf2N-] at 303 K that differ from experimental data by 25%, on average.21

7224


MAPS software and to the reviewers for very constructive suggestions that improved the quality of the paper. Supporting Information Available: Atomic numbering of ionic liquids, additional parameters and specifications, and electrostatic charges. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Wasserscheid, P.; Welton, T., Eds. Ionic Liquids in Synthesis, 2nd ed.; Wiley-VCH Verlag: Weinheim, Germany, 2007. (2) Earle, M. J.; Seddon, K. R. Pure Appl. Chem. 2000, 72, 1391. (3) Brennecke, J. F.; Maginn, E. J. AIChE J. 2001, 47, 2384. (4) Olivier-Bourbigou, H.; Magna, L. J. Mol. Catal. A: Chem. 2002, 182-183, 419. (5) Tempel, D. J.; Henderson, P. B.; Brzozowksi, J. R.; Pearlstein, R. M.; Cheng, H. S. J. Am. Chem. Soc. 2008, 130, 400. (6) Earle, M. J.; Esperanç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. Nature (London) 2006, 439, 831. (7) Cadena, C.; Anthony, J. L.; Shah, J. K.; Morrow, T. I.; Brennecke, J. F.; Maginn, E. J. J. Am. Chem. Soc. 2004, 126, 5300. (8) Huang, X.; Margulis, C. J.; Li, Y.; Berne, B. J. J. Am. Chem. Soc. 2005, 127, 17842. (9) Morrow, T. I.; Maginn, E. J. J. Phys. Chem. B 2002, 106, 12807. (10) Raabe, G.; Ko¨hler, J. J. Chem. Phys. 2008, 128, 154509. (11) Rey-Castro, C.; Tormo, A. L.; Vega, L. F. Fluid Phase Equilib. 2007, 256, 62. (12) Qiao, B.; Krekeler, C.; Berger, R.; Delle Site, L.; Holm, C. J. Phys. Chem. B 2008, 112, 1743. (13) Liu, Z.; Huang, S.; Wang, W. J. Phys. Chem. B 2004, 108, 12978. (14) Canongia Lopes, J. N.; Deschamps, J.; Pa´dua, A. A. H. J. Phys. Chem. B 2004, 108, 2038. ´ .P.S.; Tuma, D.; Maurer, G. J. Chem. (15) Kumełan, J.; Kamps, A Thermodyn. 2006, 38, 1396. (16) Tokuda, H.; Hayamizu, K.; Ishii, K.; Susan, M. A. B. H.; Watanabe, M. J. Phys. Chem. B 2005, 109, 6103. (17) Santos, L. M. N. B. F.; Canongia Lopes, J. N.; Coutinho, J. A. P.; Esperança, J. M. S. S.; Gomes, L. R.; Marrucho, I. M.; Rebelo, L. P. N. J. Am. Chem. Soc. 2007, 129, 284. (18) Kelkar, M. S.; Magginn, E. J. J. Phys. Chem. B 2007, 111, 9424. (19) Canongia Lopes, J. N.; Pa´dua, A. A. H. J. Phys. Chem. B 2004, 108, 16893. (20) Canongia Lopes, J. N.; Pa´dua, A. A. H. J. Phys. Chem. B 2006, 110, 7485. (21) Ko¨ddermann, T.; Paschek, D.; Ludwig, R. ChemPhysChem 2007, 8, 2464. (22) Ko¨ddermann, T.; Paschek, D.; Ludwig, R. ChemPhysChem 2008, 9, 549. (23) Canongia Lopes, J. N. A.; Pa´dua, A. A. H. J. Phys. Chem. B 2006, 110, 3330. (24) Wang, Y.; Voth, G. A. J. Am. Chem. Soc. 2005, 127, 12192. (25) Jiang, W.; Wang, Y.; Yan, T.; Voth, G. A. J. Phys. Chem. B 2008, 112, 1132. (26) Margules, C. J. Mol. Phys. 2004, 102, 829. (27) Urahata, S. M.; Ribeiro, M. C. C. J. Chem. Phys. 2005, 122, 024511. (28) Triolo, A.; Russina, O.; Bleif, H.-J.; Di Cola, E. J. Phys. Chem. B 2007, 111, 4641. (29) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li,

Logotheti et al. X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision B.05; Gaussian, Inc.: Wallingford, CT, 2004. (30) Becke, A. D. J. Chem. Phys. 1992, 96, 2155. (31) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (32) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361. (33) More information can be found at: www.scienomics.com. (34) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. J. Comput. Chem. 2005, 26, 1781. (35) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990. (36) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117. (37) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. J. Chem. Phys. 1995, 103, 4613. (38) Kale, L.; Skeel, R.; Bhandarkar, M.; Brunner, R.; Gursoy, A.; Krawetz, N.; Phillips, J.; Shinozaki, A.; Varadarajan, K.; Schulten, K. J. Comput. Phys. 1999, 151, 283. (39) Cadena, C.; Zhao, Q.; Snurr, R. Q.; Maginn, E. J. J. Phys. Chem. B 2006, 110, 2821. (40) Cadena, C.; Maginn, E. J. J. Phys. Chem. B 2006, 110, 18026. (41) Jeon, Y.; Sung, J.; Seo, C.; Lim, H.; Cheong, H.; Kang, M.; Moon, B.; Ouchi, Y.; Kim, D. J. Phys. Chem. B 2008, 112, 4735. (42) Potoff, J. J.; Siepmann, J. I. AIChE J. 2001, 47, 1676. (43) Caporali, S.; Bardi, U.; Lavacchi, A. J. Electron Spectrosc. Relat. Phenom. 2006, 151, 4. (44) Gomes de Azevedo, R.; Esperanca, J. M. S. S.; Szydlowski, J.; Visak, Z. P.; Pires, P. F.; Guedes, H. J. R.; Rebelo, L. P. N. J. Chem. Thermodyn. 2005, 37, 888. (45) Esperanca, J. M. S. S.; Visak, Z. P.; Plechkova, N. V.; Seddon, K. R.; Guedes, H. J. R.; Rebelo, L. P. N. J. Chem. Eng. Data 2006, 51, 2009. (46) Berg, R. W.; Deetlefs, M.; Seddon, K. R.; Shim, I.; Thompson, J. M. J. Phys. Chem. B 2005, 109, 19018. (47) de Andrade, J.; Boës, E. S.; Stassen, H. J. Phys. Chem. B 2002, 106, 3546. (48) Takahashi, S.; Suzuya, K.; Kohara, S.; Koura, N.; Curtiss, L. A.; Saboungi, M. L. Z. Phys. Chem (Munich) 1999, 209, 209. (49) Hardacre, C.; Holbrey, J. D.; McMath, S. E. J.; Bowron, D. T.; Soper, A. K. J. Chem. Phys. 2003, 118, 273. (50) Hardacre, C.; McMath, S. E. J.; Nieuwenhuyzen, M.; Bowron, D. T.; Soper, A. K. J. Phys.: Condens. Matter 2003, 15, S159. (51) Shah, J. K.; Brennecke, J. F.; Maginn, E. J. Green Chem. 2002, 4, 112. (52) Margulis, C. J.; Stern, H. A.; Berne, B. J. J. Phys. Chem. B 2002, 106, 12017. (53) Hanke, C. G.; Price, S. L.; Lynden-Bell, R. M. Mol. Phys. 2001, 99, 801. (54) Del Popolo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 1744. (55) Doxastakis, M.; Theodorou, D. N.; Fytas, G.; Kremer, F.; Faller, R.; Mu¨ller-Plathe, F.; Hadjichristidis, N. J. Chem. Phys. 2003, 119, 6883. (56) Logotheti, G. E.; Theodorou, D. N. Macromolecules 2007, 40, 2235. (57) Arrighi, V.; Batt-Coutrot, D.; Zhang, C.; Telling, M. T. F.; Triolo, A. J. Chem. Phys. 2003, 119, 1271. (58) Wang, Y.; Jiang, W.; Yan, T.; Voth, G. A. Acc. Chem. Res. 2007, 40, 1193.

JP806999Y

Molecular Modeling of Imidazolium-Based [Tf2Nâ] Ionic Liquids

Recommend Documents

Molecular Modeling of Imidazolium-Based [Tf2Nâ] Ionic Liquids