Systematic Refinement of Canongia Lopes–Pádua Force Field for

Article pubs.acs.org/JPCB

Systematic Refinement of Canongia Lopes−Pádua Force Field for Pyrrolidinium-Based Ionic Liquids Vitaly V. Chaban*,† and Iuliia V. Voroshylova‡ †

Instituto de Ciência e Tecnologia, Universidade Federal de São Paulo, 12247-014 São José dos Campos, SP, Brazil CIQ/REQUIMTE − Department of Chemistry and Biochemistry, Faculty of Sciences, University of Porto, Rua do Campo Alegre 687, 4169-007 Porto, Portugal

‡

S Supporting Information *

ABSTRACT: Reliable force field (FF) is a central issue in successful prediction of physical chemical properties via computer simulations. While Canongia LopesPád ua (CL&P) FF provides good to excellent thermodynamics and structure of pure room-temperature ionic liquids (RTILs), it suffers from drastically and systematically underestimated ionic motion. This occurs due to neglected partial electron transfer from the anion to the cation, resulting in unphysically small simulated self-diffusion and conductivity and high shear viscosities. We report a systematic refinement of the CL&P FF for six pyrrolidinium-based RTILs (1-N-butyl-1-methylpyrrolidinium dicyanamide, triflate, bis(fluorosulfonyl)imide, bis(trifluoromethanesulfonyl)imide, tetrafluoroborate, chloride). The elaborated procedure accounts for specific cation−anion interactions in the liquid phase. Once these interactions are described effectively, experimentally determined transport properties can be reproduced with an acceptable accuracy. Together with the original CL&P parameters, our force field fosters computational investigation of ionic liquids. In addition, the reported results shed more light on the chemical nature of cation−anion binding in various families of RTILs.

■

pioneered force field development for ionic liquids.3−6 A large variety of RTILs have been parametrized within the framework of their well-established methodology (CL&P).2−11 The developed models are internally consistent, transferrable, and compatible,2 providing a good to excellent accuracy of the most simulated thermodynamic and structure properties at room conditions and somewhat elevated temperatures. The major advantage of the CL&P FF over other FFs is its availability for a large number of ionic species. Another comprehensive contribution to the RTIL simulations has been delivered by Sambasivarao and Acevedo.12 Sixty-eight unique combinations of the room-temperature ionic liquids featuring 1-alkyl-3-methylimidazolium, N-alkylpyridinium, and choline cations, supplemented by an extensive set of anions, have been parametrized using harmonic bonded forces and additive pairwise potential functions. The developed models have been tested using the Metropolis Monte Carlo simulations. Reduced charges on the ions of RTILs have been implemented for the first time by Morrow and Maginn for molecular dymanics simulation of 1-n-butyl-3-methylimidazolium hexafluorophosphate.13 Ab initio calculations have yielded total charges on the cation and anion of +0.904 and −0.904,

INTRODUCTION Computer simulation has become a mighty and widely recognized tool to predict physical chemical properties of room-temperature ionic liquids (RTILs).1,2 A basic search performed in the ISI Web of Knowledge database on November 7, 2014, reveals that 373 papers were devoted to computer simulations of RTILs only during 2014. Compare this number with a total number of papers devoted to ionic liquids, 7072, during 2014. Hereby, 5.3% constitute a very significant portion, provided that an average computer simulation contribution is more intellectually expensive than a direct physical chemical experiment employing an established methodology. Most numerical simulations report atomisticresolution molecular dynamics, based on pairwise preparametrized interaction potential functions. A smaller number of works employ more sophisticated interaction potentials, Monte Carlo techniques, zero-temperature electronic structure studies, ab initio molecular dynamics simulations, coarse-grained representations of interacting sites, etc. Reliable sampling of phase space and correct representation of interatomic binding functions (known as force field, FF) are a central issue in successful prediction of physical chemical properties via various sorts of computer simulation techniques. Several developments for the state-of-the-art atomistic description of ionic liquids have been introduced over the recent years. Canongia Lopes and Pádua have essentially © 2015 American Chemical Society

Received: December 9, 2014 Revised: March 29, 2015 Published: March 31, 2015 6242

DOI: 10.1021/jp5122678 J. Phys. Chem. B 2015, 119, 6242−6249

Article

The Journal of Physical Chemistry B respectively. The results of molecular dynamics simulation has been obtained in good (for infrared spectra, volume expansivity, and isothermal compressibility) to excellent (for molar volumes) agreement with experimental data. In turn, Yan et al.14 have applied Thole model to produce polarizable ions for ionic liquid 1-ethyl-3-methylimidazolium nitrate. Shear viscosity and diffusion constants calculated at 400 K using both polarizable and nonpolarizable versions of the model have shown that the effects of polarizability are quite substantial. The polarizable model results are in better agreement with the experimental values. The discussion has been extended on the enhanced dynamical properties later on.15 It has been shown that the polarizable model consistently exhibits faster relaxations, higher self-diffusion and ionic conductivity, and lower shear viscosity than the nonpolarizable model. Recently, Borodin has proposed a polarizable model for a significant number of RTILs.16 A set of simulated thermodynamic (density, heat of vaporization) and transport properties (selfdiffusion coefficients, shear viscosity, ionic conductivity) appear in a reasonable agreement16 with the available experimental results. We have previously made certain contribution to the FF development and refinement by coupling electronic structure and molecular mechanical descriptions of interionic interactions.17−21 A united-atom FF for pyrrolidinium RTILs with bis(trifluoromethanesulfonyl)imide and bis(fluorosulfonyl)imide anions22,23 has been reported recently. The developed model has been used to carry out molecular dynamics simulation studies of the structure and the differential capacitance near graphite electrode. However, all-atom FF provides physical insights with higher precision. The present work further elaborates/justifies an approach to develop efficient atomistic force field models for RTILs under a simple approximation of additive nonbonded interactions. Starting from the FF for isolated ions in vacuum (CL&P), we consider two major corrections: (1) for electron transfer effects between the cation and the anion and (2) for cation− anion atom−atom closest-approach distances. Both corrections are driven by the results of hybrid density functional theory (DFT) for ion clusters. We apply our methodology to six RTILs representing the pyrrolidinium family: (1) 1-N-butyl-1methylpyrrolidinium tetrafluoroborate [PY14][BF4]; (2) 1-Nbutyl-1-methylpyrrolidinium bis(trifluoromethanesulfonyl)imide [PY14][TFSI]; (3) 1-N-butyl-1-methylpyrrolidinium dicyanamide [PY14][DCA]; (4) 1-N-butyl-1-methylpyrrolidinium bis(fluorosulfonyl)imide [PY14][FSI]; (5) 1-N-butyl-1methylpyrrolidinium triflate [PY14][TF]; (6) 1-N-butyl-1methylpyrrolidinium chloride [PY14][Cl]. Figure 1 introduces ion aggregates (two ion pairs), which constitute primary building blocks of the investigated RTILs. Certain experimental data for four of these RTILs are available (see Table S1).24−43 These data can be used to evaluate the accuracy of the derived FF. We have not found any experimental properties for [PY14][BF4] and [PY14][Cl]. Experimental efforts to address transport properties of the two mentioned RTILs are very welcome. Since 1-ethyl-1-methylpyrrolidinium tetrafluoroborate is a plastic crystal at ambient conditions,44 it is quite likely that [PY14][BF4] exhibits a higher boiling point than other pyrrolidinium-based ILs. Thus, transport properties reported for [PY14][BF4] must be taken with an appropriate caution.

Figure 1. Optimized geometries for the isolated clusters composed of four RTIL ions each: two 1-ethyl-1-methylpyrrolidinium cations and two respective anions. 1-Ethyl-1-methylpyrrolidinium cations are used in the ab initio calculations instead of 1-N-butyl-1-methylpyrrolidinium cation.

electron density using the recently proposed high-quality hybrid exchange-correlation functional, omega B97X-D.45 The exchange energy is combined with the exact energy from the Hartree−Fock theory. Empirical atom−atom dispersion correction is included. According to Chai and Head-Gordon, the omega B97X-D functional simultaneously yields an improved accuracy for thermochemistry, electron kinetics, and noncovalent interactions.45 We construct the wave function using a comprehensive correlation-consistent basis set proposed by Dunning and co-workers.46 Based on our preliminary tests, triple-ζ basis set from this family, augmented with the diffuse functions (aug-cc-pVTZ),46 provides an accurate description of the electronic density and energy levels. Certain known defects, such as excessive valence electron delocalization, widely observed in the case of pure DFT methods, were not reported for the presently selected method and basis set. Furthermore, we investigate how electron transfer between the anion and the cation is impacted by the number of ion pairs. We test ion clusters containing 1, 2, 3, and 4 ions pairs (i.e., up to 8 ions). Since the resulting clusters are essentially large in terms of electrons, we employed a smaller basis set, 6-31G(d), for this particular benchmark (Figure 2). Electrostatic potential (ESP), derived from the electronic wave function, was fitted using a set of point charges localized on each atom of the cation and anion, including hydrogen atoms. In the case of ion pairs (i.e., neutral systems), an additional constraint was imposed during the fitting procedure to rigorously reproduce dipole moment of the ionic aggregate. Atomic radii were assigned on the basis of the CHELPG scheme.47 All electronic structure calculations and geometry optimizations of ionic clusters were performed using the GAUSSIAN 09, revision D, program. The thermodynamics (heat of vaporization, ΔHvap, mass density, d) and transport properties (self-diffusion coefficients, D±, shear viscosity, η) have been obtained using equilibrium and nonequilibrium molecular dynamics (MD) simulations, as detailed below. The basic list of the simulated systems is

■

METHODS AND TOOLS Cation−anion specific interactions were investigated using the hybrid electronic density theory (Figure 1). We describe an 6243

DOI: 10.1021/jp5122678 J. Phys. Chem. B 2015, 119, 6242−6249

Article

The Journal of Physical Chemistry B

calculation. The standard error of each viscosity value was evaluated from statistical processing of 0−15, 15−30, and 30− 45 ns trajectory fractions. Hereby, these trajectory fractions were treated as independent measurements in statistical analysis. All systems were simulated in the constant-pressure constanttemperature ensemble. The equations-of-motion were propagated with a time step of 1.0 fs. The electrostatic interactions were simulated using direct Coulomb law up to 1.3 nm of separation between the interaction sites. The electrostatic interactions beyond 1.3 nm were accounted for by computationally efficient particle-mesh Ewald (PME) method.49 It is important to use the PME method in the case of ionic systems, since electrostatic energy beyond the cutoff usually contributes 40−60% of total electrostatic energy. We used a smaller realspace cutoff, as compared to the original works of Canongia Lopes−Pádua.5,6,8,9 This methodological adjustment does not change the accuracy noticeably, according to our preliminary assessment, but allows to speed up the MD calculations. The Lennard-Jones 12−6 interactions were smoothly brought down to 0 from 1.1 to 1.2 nm using the classical shifted force technique. The constant temperature (298, 313, 338, and 353 K) was maintained by the Bussi−Donadio−Parrinello velocity rescaling thermostat (with a time constant of 0.5 ps),50 which provides a correct velocity distribution for a statistical mechanical ensemble. The constant pressure was maintained by Parrinello−Rahman barostat51 with a time constant of 4.0 ps and a compressibility constant of 4.5 × 10−5 bar−1. All molecular dynamics trajectories were propagated using the GROMACS simulation suite.49,52,53 Analysis of properties was done using the in-home tools and supplementary utilities distributed with the GROMACS package.49,52,53

Figure 2. Optimized ionic configurations of [PY12]n[DCA]n for n = 1, 2, 3, and 4.

provided in Table 2. The ions of RTILs were placed in the cubic periodic MD boxes, whose densities were calculated to correspond to ambient pressure at the requested temperature (298, 313, 338, and 353 K). The simulated temperatures were selected in view of the currently available experimental data and a primary interest for existing and prospective applications. The first 10 ns of recorded trajectory for every system was regarded as equilibration, whereas subsequent trajectory part was used to sample the phase space and derive physical chemical properties of interest (see Table 2 for simulation durations). Note that each system at each temperature underwent nonequilibrium simulation (for shear viscosity) and equilibrium simulation (for all other properties). Cartesian coordinates were saved every 1.000 ps, and data on energy dissipation (for viscosity calculation) were saved every 0.01 ps. More frequent saving of trajectory components was preliminarily tested, but no systematic accuracy improvement was found. Self-diffusion coefficients were computed from meansquare displacements of atomic positions versus time. Prior to these calculations, the trajectories must be preprocessed to remove information about periodic boundary conditions. Shear viscosity was calculated using cosine-shape acceleration of all ions simultaneously. The mathematical foundation of this method and its applications to molecular dynamics simulations were, in great detail, discussed by Hess.48 The first 5 ns of the simulation was used for an accelerated ionic flow to be established. The subsequent 45 ns was used for the viscosity

■

FORCE FIELD DERIVATION 1-Ethyl-1-methylpyrrolidinium cation was used instead of 1-Nbutyl-1-methylpyrrolidinium cation for electronic structure calculations. The reason for such a solution is that electrostatically neutral and nonpolar −CH2−CH3 fragment does not participate in charge transfer between the anion and the cation. Avoiding −CH2−CH3 significantly decreases computational cost of the system and hastens convergence of self-consistent field procedures. Table 1 lists point charges obtained through fitting of electrostatic potential at the surface of the isolated 1-ethyl-1methylpyrrolidinium cation and six neutral ion pairs containing 1-ethyl-1-methylpyrrolidinium cation and various anions. The

Table 1. Electrostatic Point Charges Condensed on the Non-Hydrogen Atoms of N-1-Ethyl-1-methylpyrrolidinium Cationa 1-ethyl-1-methylpyrrolidinium+ + anion −

atom

isolated

[BF4]

[Cl]−

[DCA]−

[FSI]−

[TF]−

[TFSI]−

C-1 (ring) C-2 (ring) C-3 (ring) C-4 (ring) N-5 C-6 (CH3) C-7 (CH2) C-8 (CH3) total

+0.19 +0.04 +0.15 +0.10 +0.12 +0.18 +0.14 +0.08 +1.00

+0.20 0.00 +0.09 +0.12 +0.13 +0.07 +0.24 −0.01 +0.84

+0.18 −0.01 +0.08 +0.10 +0.19 +0.05 +0.19 −0.01 +0.77

+0.12 +0.05 +0.11 +0.06 +0.08 +0.09 +0.28 −0.02 +0.77

+0.20 −0.01 +0.10 +0.09 +0.18 +0.08 +0.23 +0.01 +0.88

+0.20 −0.01 +0.09 +0.11 +0.13 +0.08 +0.25 0.00 +0.88

+0.24 −0.06 +0.14 +0.07 +0.09 +0.11 +0.28 −0.01 +0.86

a

The charges were assigned according to the CHELPG scheme. An isolated cation state is compared to ion pairs containing various anions. The electronic structure calculations for all systems were conducted at the omega B97XD/aug-cc-pVTZ level of theory.45,46 6244

DOI: 10.1021/jp5122678 J. Phys. Chem. B 2015, 119, 6242−6249

Article

The Journal of Physical Chemistry B Table 2. List of Systems Simulated in the Present Worka sampling times, ns no.

RTIL

no. of ion pairs

no. of interaction centers

equilibriumb

nonequilibriumc

1−4 5−8 9−12 13−16 17−20 21−24

[PY14][TFSI] [PY14][BF4] [PY14][DCA] [PY14][FSI] [PY14][TF] [PY14][Cl]

125 150 150 125 125 175

5625 5250 5250 4875 4750 5425

20 30 20 25 20 30

50 70 50 50 50 70

a

Each ionic composition was simulated at 298, 313, 338, and 353 K (i.e., four MD systems per ionic composition). The quantity of ion pairs per system was selected with respect to the cation and the anion sizes. The durations of trajectories for various RTILs were selected in a way so that computational cost per system was roughly equal in all cases. Note that required sampling strongly depends on the system temperature and, ultimately, on the average self-diffusion constant in the respective system. bEquilibrium MD simulation was used to sample mass density, heat of vaporization, and self-diffusion coefficients of the cation and the anion. cNonequilibrium MD simulation was used for shear viscosity calculation as described above.

clusters with [DCA]− and [Cl]− as a function of the cluster size. Addition of the second ion pair brings insignificant decrease of the total charge, ca. 1%, while addition of more ions does not exhibit a systematic impact. It should be noted that electrostatic potential is well-defined only at the surface of the ion cluster; therefore, further addition of ions does not necessarily correspond to a more realistic description of point charges. In turn, implicit solvent models may capture polarization effect of the surrounding solvent but does not describe an electron transfer phenomenon (since the solvent is implicit). It is questionable whether implementation of periodic boundary conditions (PBC) is more beneficial in the present case, since PBC does not represent a long-range structure of RTIL. Because of small size of the unit cell, reproduction of experimental densitywhether in the constant-volume ensemble or in the constant-pressure ensembleappears problematic. A few works have been reported recently address electronic polarization in RTILs as a mean field54 with optical dielectric contant.55 The scaling factors proposed in this work are between ±0.7−0.8 per ion. These values are similar to 1/√2. Consequently, optical dielectric constant due to electronic polarization equals to ca. 2, which is in concordance with refs 54 and 55. Such a similarity adds more physical background to our procedure. Since the total charge of an ion does not change significantly in larger clusters, both the result for one ion pair and two ion pairs can be used. In this work, we chose the result for two ion pairs (four ions), averaged to two decimal digits, to be used for the force field derivation. It should be noted that various DFT functionals provide reasonably good mutual agreement. Compare, in the case of [PY12]2[DCA]2, PBE56 provides ±0.69e/ion, B97D45 provides ±0.70e/ion, M06L57 provides ±0.72e/ion, and BLYP58provides ±0.70e/ion. Although the difference is somewhat larger than the difference between various number of ion pairs, it can be considered sufficiently small in the context of force field derivation. We used the charges from the PBE functional56 for the force field derivation. The present choice is arbitrary. The set of immediate charges, such as those summarized in Table 1, cannot be directly transferred to the classical force field. All of these point charges are small. They fluctuate drastically upon thermal motion of RTIL. Instead, the symmetry of charge distribution in CL&P FF can be preserved, whereas an electrostatic potential can be decreased through uniform scaling of all charges of the pyrrolidine ring of the cation and the anion. This solution was successfully applied in

point charges on the hydrogen atoms were summed up with the charge on the carbon atom for simplicity of representation. All ion pairs exhibit a significant electron delocalization, indicatingamong other observationsstrong cation−anion binding in these configurations. The total charges of the respective anions are the same as of 1-ethyl-1-methylpyrrolidinium, but with an opposite sign. Cation−anion binding breaks symmetries of both isolated cation and anion. The point charges on the same heavy atoms of the cation (carbon, nitrogen) differ very significantly when the cation is paired with different anions. Compare, q(N) = +0.12e in the isolated [PY12]+, q(N) = +0.19e in [PY12][Cl], and q(N) = +0.08e in [PY12][DCA]. Note that these calculations do not account for thermal motion effects. Thermal motion presumably brings even larger fluctuations of point charges. Since the charge of C8 (CH3 radical, which belongs to the ethyl chain, Table 1) appears nearly zero with all anions, we disregard it in the further calculations of larger ion clusters (Figure 1). The total charge on the ion can, in principle, depend on the size of ion cluster. It is a probable situation in compounds, where more than one coordination center per molecule (ion) is present. In such cases, a single counterion does not occupy all possible sites, which are realized in the condensed phase. Hence, charge transfer may be described partially but not completely. Figure 2 depicts optimized geometries of [PY12][DCA] clusters containing 2−8 RTIL ions (1−4 ion pairs). Figure 3 summarizes total charge on the [PY12] cation in its

Figure 3. Electrostatic (CHELPG) charge localized on the [PY12]+ cation in the [PY12]n[DCA]n and [PY12]n[Cl]n, n = 1, 2, 3, and 4. The depicted charge was computed as a sum of all point electrostatic charges on the [PY12]+ atoms, including hydrogen atoms. The case of zero ion pairs corresponds to an isolated cation. Its charge equals to unity according to the number of electrons in the considered nuclei− electron system. 6245

DOI: 10.1021/jp5122678 J. Phys. Chem. B 2015, 119, 6242−6249

Article

The Journal of Physical Chemistry B our previous works. The so-called scaling factor equals to the total charge on an ion suggested by DFT. Remember that electrostatically neutral moiety of the hydrocarbon chain, −CH2−CH3, does not participate in charge transfer, and hence, charge scaling is not necessary for it. We implement uniform scaling factors for [PY14][DCA], [PY14][FSI], [PY14][TF], [PY14][BF4], and [PY14][Cl]. It could have been also applied with reasonable success to [PY14][TFSI], as illustrated before.17 However, in view of the two sulfur atoms in [TFSI]−, we attempted to account for their different polarizability by transferring charges from an ab initio calculation of [PY12]2[TFSI]2. Direct transfer of the ESP point charges from the ab initio calculation to the Coulomb potential is not doable because formation of an ion pair breaks an original symmetry of isolated ions. That is, symmetric atoms of sulfur, oxygen, carbon, and fluorine in the [TFSI]− anion obtain significantly alternating charges. Being averaged with respect to the first decimal digit, the point charges are −0.4e (N), +0.6e (S), −0.4e (O), +0.3e (C), and −0.1e (F), whereas the charges of the 1butyl-1-methylpyrrolidinium cation were uniformly scaled. It should be noted that we did not tune the bonded parameters in the CL&P FF, since predicted conformations are considered correct. We performed a few tests to ensure that decreased partial charges do not induce drastic changes in the cationic and anionic conformations All proposed charges for all six ionic liquid are summarized in the FF files provided in the Supporting Information. Unlike in the pyridinium-based RTILs, the cation−anion closest-approach distance in the pyrrolidinium-based RTILs is described well by the existing CL&P FF parameters. Anions are coordinated by the methyl group, as opposed to coordination by the ring. Coordination by ring occurs in the pyridinium-based RTILs. This is, in our opinion, a principal difference between pyridinium-based and pyrrolidinium-based RTILs, which was not underlined in the literature thus far. The refined FF treats all atoms explicitly, including hydrogen atoms. This is done for compatibility with the CL&P FF.2,5,6,8−11 Construction of the united atom force field is an interesting option, which would allow avoiding a number of dihedral functions in the Hamiltonian and therefore enhance sampling.

Figure 4. Heat of vaporization for the six simulated pyrrolidiniumbased ionic liquids at 298, 313, 338, and 353 K.

■

FORCE FIELD VALIDATION Figure 4 depicts temperature dependence of heat of vaporization, Hvap. All curves are nearly linear over the 298−353 K temperature range. In the case of molecular liquids, heat of vaporization is normally measured at the normal boiling point and then extrapolated to 298 K. The difference obtained after extrapolation is often inferior to the uncertainty of the measured value, which is high. The decrease of Hvap with temperature increase is modest over the condensed matter temperature range. However, Hvap decays drastically as temperature approaches critical point, which is not the case in the present study. The experimental Hvap of [PY14][DCA] amounts to 161 kJ mol−1 (Table S1), which is in satisfactory agreement with the simulated Hvap = 139 kJ mol−1. The simulated Hvap of [PY14][TFSI] is equal to 150 kJ mol−1, whereas the experimental values are 152 and 144.4 kJ mol−1 (Table S1). Thus, the developed model performs perfectly. Mass densities of the six pyrrolidinium-based RTILs (Figure 5) are in concordance with the experimental results (Table S1).24−43 Compare, 958 vs 950 kg m−3 for [PY14][DCA], 1326 vs 1302 kg m−3 for [PY14][FSI], and 1398 vs 1399 kg m−3 for

Figure 5. Mass density for the six simulated pyrrolidinium-based ionic liquids at 298, 313, 338, and 353 K. The size of circles was selected to be larger than the corresponding block-averaged standard errors (from simulations).

[PY14][TFSI]. The case of [PY14][TF] is the least successful out of these examples, 1204 vs 1256 kg m−3. The discrepancy of 4% is within acceptable boundaries for classical molecular dynamics simulations. The reason for this discrepancy may be associated with less accurate atomic radii for specific cation− anion binding. [PY14][DCA] with a low density of 958 kg m−3 at 298 K constitutes an interesting case. Less dense substances, in general, feature higher ionic conductivity and smaller shear 6246

DOI: 10.1021/jp5122678 J. Phys. Chem. B 2015, 119, 6242−6249

Article

The Journal of Physical Chemistry B viscosity. In turn, [PY14][TFSI] represents a notably dense ionic liquid, thanks to relatively heavy sulfur and fluorine atoms in the anion. Indeed, both cation and anion of [PY14][DCA] exhibit higher self-diffusivities, 19 × 10−11 and 25 × 10−11 m2 s−1 (Figure 6) than any other ion from our list. These diffusions

Figure 7. Shear viscosity for the six simulated pyrrolidinium-based ionic liquids at 298, 313, 338, and 353 K. Figure 6. Self-diffusion constants for the six simulated pyrrolidinium ionic liquids at 298, 313, 338, and 353 K. The self-diffusion constants of the [PY14]+ cation are depicted in the upper panel, whereas the self-diffusion constants of the respective anions are depicted in the lower panel. The size of symbols is chosen to exceed standard errors of the computed self-diffusion constants at each temperature.

constants are in excellent agreement with the available experiments (Table S1). The slowest diffusion is observed in [PY14][FSI] and [PY14][TFSI], which is in line with their molecular masses and densities. In no RTIL, self-diffusion of cation differs from that of counterion significantly. Therefore, the case of very mobile cation and relatively slow anion or contrariwise is not observed. Small anions (chloride, tetrafluoroborate) coupled with [PY14]+ become immobile due to a strong electrostatic attraction (high potential energies) in these systems. The simulated shear viscosities at 298, 313, 338, and 353 K are depicted in Figure 7. The comparison between experimental24−43 and simulated viscosities is provided in Figure 8. Exponential growth of viscosity with linear temperature increase was successfully achieved. The highest viscosities were observed for [PY14][FSI] (118 cP) and [PY14][BF4] (375 cP) at room temperature. In the first case, it predominantly occurs due to a large anionic species, whereas strong electrostatic contribution (polar covalent boron− fluorine bonds) is responsible for the second case. Interestingly, the chloride liquid is not most viscous in the current data set. According to the developed force field and DFT calculations, the interaction between [Cl]− and [PY14]+ is not extremely strong. Indeed, Figure 1 suggests that [Cl]− is coordinated by the hydrocarbon (methyl) groups of the two neighboring pyrrolidinium-based cations, as opposed to coordination by the pyrrolidine ring. Unfortunately, experimental data are not yet available for the chloride containing RTIL. Correspondence between simulated and experimental data24−43 is best at 353 K (Figure 8). We have ensured that poorer performance of the developed FF at lower temperatures (298−338 K) is not due to insufficient sampling. This test has

Figure 8. Divergence of the simulated shear viscosities and the experimental data.24−31,33−43,59 The divergence (in %) is obtained as η = [(η_sim − η_exp)/η_exp] × 100%.

been done by computing viscosities from the subsequent parts of the trajectory and comparing them to one another. No drift was observed. The computed viscosities of [PY14][FSI] at all temperatures appear significantly and systematically higher than the experimental values (Table S1). The physical/computational reasons of such a behavior are not clear. It is especially surprising because other models (five pure liquids) of pyrrolidinium-based RTILs exhibit an acceptable performance. Our FF derivation philosophy does not allow changing any interaction parameters without a clear physical implication. Therefore, we retain our FF parameters as they are. One must not exclude a possibility that the experimental viscosities are underestimated. Such a situation can occur, for instance, due to poor drying of this particular RTIL prior to viscosity measurements. The reported shear viscosities were computed using separate nonequilibrium MD simulations based on how fast extra energy is dissipated in the system. The deviation from equilibrium 6247

DOI: 10.1021/jp5122678 J. Phys. Chem. B 2015, 119, 6242−6249

Article

The Journal of Physical Chemistry B

(5) Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. A. H. Modeling Ionic Liquids Using a Systematic All-Atom Force Field. J. Phys. Chem. B 2004, 108 (6), 2038−2047. (6) Canongia Lopes, J. N.; Pádua, A. A. H. Molecular Force Field for Ionic Liquids Composed of Triflate or Bistriflylimide Anions. J. Phys. Chem. B 2004, 108 (43), 16893−16898. (7) Canongia Lopes, J. N. A.; Pádua, A. A. H. Using Spectroscopic Data on Imidazolium Cation Conformations To Test a Molecular Force Field for Ionic Liquids. J. Phys. Chem. B 2006, 110 (14), 7485− 7489. (8) Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. A. H. Modeling Ionic Liquids of the 1-Alkyl-3-methylimidazolium Family Using an AllAtom Force Field. Ionic Liquids IIIa: Fundamentals, Progress, Challenges, and Opportunities, Properties and Structure 2005, 901, 134−149. (9) Canongia Lopes, J. N.; Pádua, 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 (39), 19586−19592. (10) Canongia Lopes, J. N.; Pádua, A. A. H.; Shimizu, K. Molecular Force Field for Ionic Liquids IV: Trialkylimidazolium and Alkoxycarbonyl-Imidazolium Cations; Alkylsulfonate and Alkylsulfate Anions. J. Phys. Chem. B 2008, 112 (16), 5039−5046. (11) Shimizu, K.; Almantariotis, D.; Gomes, M. F. C.; Pádua, A. A. H.; Canongia Lopes, J. N. Molecular Force Field for Ionic Liquids V: Hydroxyethylimidazolium, Dimethoxy-2-Methylimidazolium, and Fluoroalkylimidazolium Cations and Bis(Fluorosulfonyl)Amide, Perfluoroalkanesulfonylamide, and Fluoroalkylfluorophosphate Anions. J. Phys. Chem. B 2010, 114 (10), 3592−3600. (12) Sambasivarao, S. V.; Acevedo, O. Development of OPLS-AA Force Field Parameters for 68 Unique Ionic Liquids. J. Chem. Theory Comput. 2009, 5 (4), 1038−1050. (13) Morrow, T. I.; Maginn, E. J. Molecular Dynamics Study of the Ionic Liquid 1-n-Butyl-3-methylimidazolium Hexafluorophosphate. J. Phys. Chem. B 2002, 106 (49), 12807−12813. (14) Yan, T.; Burnham, C. J.; Del Pópolo, M. G.; Voth, G. A. Molecular Dynamics Simulation of Ionic Liquids: The Effect of Electronic Polarizability. J. Phys. Chem. B 2004, 108 (32), 11877− 11881. (15) Yan, T.; Wang, Y.; Knox, C. On the Dynamics of Ionic Liquids: Comparisons between Electronically Polarizable and Nonpolarizable Models II. J. Phys. Chem. B 2010, 114 (20), 6886−6904. (16) Borodin, O. Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids. J. Phys. Chem. B 2009, 113 (33), 11463−11478. (17) Chaban, V. Polarizability versus Mobility: Atomistic Force Field for Ionic Liquids. Phys. Chem. Chem. Phys. 2011, 13 (35), 16055− 16062. (18) Chaban, V.; Maciel, C.; Fileti, E. Does the Like Dissolves Like Rule Hold for Fullerene and Ionic Liquids? J. Solution Chem. 2014, 43 (6), 1019−1031. (19) Chaban, V. V.; Voroshylova, 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 (17), 7910−7920. (20) Chaban, V. V.; Voroshylova, I. V.; Kalugin, O. N.; Prezhdo, O. V. Acetonitrile Boosts Conductivity of Imidazolium Ionic Liquids. J. Phys. Chem. B 2012, 116 (26), 7719−7727. (21) Voroshylova, I. V.; Chaban, V. V. Atomistic Force Field for Pyridinium-Based Ionic Liquids: Reliable Transport Properties. J. Phys. Chem. B 2014, 118 (36), 10716−10724. (22) Vatamanu, J.; Borodin, O.; Smith, G. D. Molecular Insights into the Potential and Temperature Dependences of the Differential Capacitance of a Room-Temperature Ionic Liquid at Graphite Electrodes. J. Am. Chem. Soc. 2010, 132 (42), 14825−14833. (23) Vatamanu, J.; Borodin, O.; Smith, G. D. Molecular Simulations of the Electric Double Layer Structure, Differential Capacitance, and Charging Kinetics for N-Methyl-N-propylpyrrolidinium Bis(fluorosulfonyl)imide at Graphite Electrodes. J. Phys. Chem. B 2011, 115 (12), 3073−3084.

(degree of nonequilibrium) must be strictly controlled. Too small accelerations would result in the unsatisfactorily slow sampling, whereas too large accelerations would lead to viscosity underestimation. After preliminary tests (see Supporting Information) the cosine acceleration value of 0.02 nm ps−2 was chosen.

■

CONCLUSIONS Starting from the CL&P force field for the six pyrrolidiniumbased ionic liquids, we applied a systematic refinement for more realistic reproduction of transport properties (cationic and anionic self-diffusion coefficients, shear viscosity). The improvement was achieved thanks to description of electrostatic potential of ion aggregates, as opposed to isolated ions in vacuum (Canongia Lopes−Pádua force field). All properties of interest are in reasonable concordance with the available experimental values. Remarkably, the coincidence with the experimental properties is better in the case of anions, which do not contain sulfur. The likely reason is that sulfur, belonging to the third period of the periodic table, exhibits significantly different polarizability, as compared to all other elements constituting these RTILs. The developed force field fosters computational investigation of room-temperature ionic liquids and promotes understanding of physical and chemical atomistic-precision interactions in the condensed state of these complicated systems.

■

ASSOCIATED CONTENT

S Supporting Information *

All force field parameters derived and tested in the manuscript; Table S1: summary of all currently available experimental data on heat of vaporization, mass density, cationic and anionic selfdiffusion constants, and shear viscosity at 298, 313, 338, and 353 K. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (V.V.C.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

V.V.C. acknowledges a research grant from CAPES under “Science Without Borders” program. Iu.V.V. acknowledges a research grant from FCT with the reference SFRH/BPD/ 97918/2013.

(1) Chaban, V. V.; Prezhdo, O. V. Ionic and Molecular Liquids: Working Together for Robust Engineering. J. Phys. Chem. Lett. 2013, 4, 1423−1431. (2) Canongia Lopes, J. N.; Pádua, A. A. H. CL&P: A Generic and Systematic Force Field for Ionic Liquids Modeling. Theor. Chem. Acc. 2012, 131 (3), 1129. (3) Canongia Lopes, J. N.; Pádua, A. A. H.; Deschamps, J. Modeling Ionic Liquids Using an Extended OPLS-AA Force Field. Abstr. Pap. Am. Chem. Soc. 2003, 226 (U-618). (4) Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. A. H. Modeling Ionic Liquids Using a Systematic All-Atom Force Field. J. Phys. Chem. B 2004, 108 (30), 11250−11250. 6248

DOI: 10.1021/jp5122678 J. Phys. Chem. B 2015, 119, 6242−6249

Article

The Journal of Physical Chemistry B

Spectra of Pyrrolidinium-Based Ionic Liquids. J. Phys. Chem. B 2012, 116 (24), 7322−7327. (41) Zaitsau, D. H.; Yermalayeu, A. V.; Emel’yanenko, V. N.; Heintz, A.; Verevkin, S. P.; Schick, C.; Berdzinski, S.; Strehmel, V. Structure− Property Relationships in ILs: Vaporization Enthalpies of Pyrrolidinium Based Ionic Liquids. J. Mol. Liq. 2014, 192 (0), 171−176. (42) Zhou, Q.; Henderson, W. A.; Appetecchi, G. B.; Montanino, M.; Passerini, S. Physical and Electrochemical Properties of N-Alkyl-Nmethylpyrrolidinium Bis(fluorosulfonyl)imide Ionic Liquids: PY13FSI and PY14FSI. J. Phys. Chem. B 2008, 112 (43), 13577−13580. (43) Zhou, Z.-B.; Matsumoto, H.; Tatsumi, K. Cyclic Quaternary Ammonium Ionic Liquids with Perfluoroalkyltrifluoroborates: Synthesis, Characterization, and Properties. Chem.Eur. J. 2006, 12 (8), 2196−2212. (44) Chen, F.; Zhu, H.; Forsyth, M. Modelling Ion-Pair Geometries and Dynamics in a 1-Ethyl-1-methylpyrrolidinium-Based Ion-Conductive Crystal. ChemPhysChem 2014, 15 (16), 3530−3535. (45) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10 (44), 6615−6620. (46) Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90 (2), 1007−1023. (47) Breneman, C. M.; Wiberg, K. B. Determining Atom-Centered Monopoles from Molecular Electrostatic Potentials. The Need for High Sampling Density in Formamide Conformational Analysis. J. Comput. Chem. 1990, 11 (3), 361−373. (48) Hess, B. Determining the Shear Viscosity of Model Liquids from Molecular Dynamics Simulations. J. Chem. Phys. 2002, 116 (1), 209− 217. (49) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: A Package for Molecular Simulation and Trajectory Analysis. Mol. Model. Annu. 2001, 7 (8), 306−317. (50) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126 (1), 014101. (51) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52 (12), 7182−7190. (52) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4 (3), 435−447. (53) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26 (16), 1701−1718. (54) Leontyev, I.; Stuchebrukhov, A. Accounting for Electronic Polarization in Non-Polarizable Force Fields. Phys. Chem. Chem. Phys. 2011, 13 (7), 2613−2626. (55) Fedorov, M. V.; Kornyshev, A. A. Ionic Liquid Near a Charged Wall: Structure and Capacitance of Electrical Double Layer. J. Phys. Chem. B 2008, 112 (38), 11868−11872. (56) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77 (18), 3865− 3868. (57) Zhao, Y.; Truhlar, D. G. A New Local Density Functional for Main-Group Thermochemistry, Transition Metal Bonding, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Phys. 2006, 125 (19), 194101. (58) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37 (2), 785−789. (59) Hayamizu, K.; Tsuzuki, S.; Seki, S.; Fujii, K.; Suenaga, M.; Umebayashi, Y. Studies on the Translational and Rotational Motions of Ionic Liquids Composed of N-Methyl-N-propylpyrrolidinium (P13) Cation and Bis(trifluoromethanesulfonyl)amide and Bis(fluorosulfonyl)amide Anions and Their Binary Systems Including Lithium Salts. J. Chem. Phys. 2010, 133, (19).

(24) Bayley, P. M.; Best, A. S.; MacFarlane, D. R.; Forsyth, M. The Effect of Coordinating and Non-Coordinating Additives on the Transport Properties in Ionic Liquid Electrolytes for Lithium Batteries. Phys. Chem. Chem. Phys. 2011, 13 (10), 4632−4640. (25) Bayley, P. M.; Lane, G. H.; Lyons, L. J.; MacFarlane, D. R.; Forsyth, M. Undoing Lithium Ion Association in Ionic Liquids through the Complexation by Oligoethers. J. Phys. Chem. C 2010, 114 (48), 20569−20576. (26) Bayley, P. M.; Lane, G. H.; Rocher, N. M.; Clare, B. R.; Best, A. S.; MacFarlane, D. R.; Forsyth, M. Transport Properties of Ionic Liquid Electrolytes with Organic Diluents. Phys. Chem. Chem. Phys. 2009, 11 (33), 7202−7208. (27) Castiglione, F.; Moreno, M.; Raos, G.; Famulari, A.; Mele, A.; Appetecchi, G. B.; Passerini, S. Structural Organization and Transport Properties of Novel Pyrrolidinium-Based Ionic Liquids with Perfluoroalkyl Sulfonylimide Anions. J. Phys. Chem. B 2009, 113 (31), 10750−10759. (28) Deyko, A.; Hessey, S. G.; Licence, P.; Chernikova, E. A.; Krasovskiy, V. G.; Kustov, L. M.; Jones, R. G. The Enthalpies of Vaporisation of Ionic Liquids: New Measurements and Predictions. Phys. Chem. Chem. Phys. 2012, 14 (9), 3181−3193. (29) Deyko, A.; Lovelock, K. R. J.; Corfield, J. A.; Taylor, A. W.; Gooden, P. N.; Villar-Garcia, I. J.; Licence, P.; Jones, R. G.; Krasovskiy, V. G.; Chernikova, E. A.; Kustov, L. M. Measuring and predicting ΔHvap(298) Values of Ionic Liquids. Phys. Chem. Chem. Phys. 2009, 11 (38), 8544−8555. (30) Emel’yanenko, V. N.; Verevkin, S. P.; Heintz, A.; Corfield, J.-A.; Deyko, A.; Lovelock, K. R. J.; Licence, P.; Jones, R. G. PyrrolidiniumBased Ionic Liquids. 1-Butyl-1-methyl Pyrrolidinium Dicyanoamide: Thermochemical Measurement, Mass Spectrometry, and ab Initio Calculations. J. Phys. Chem. B 2008, 112 (37), 11734−11742. (31) Frömling, T.; Kunze, M.; Schönhoff, M.; Sundermeyer, J.; Roling, B. Enhanced Lithium Transference Numbers in Ionic Liquid Electrolytes. J. Phys. Chem. B 2008, 112 (41), 12985−12990. (32) Hayamizu, K.; Tsuzuki, S.; Seki, S.; Fujii, K.; Suenaga, M.; Umebayashi, Y. Studies on the Translational and Rotational Motions of Ionic Liquids Composed of N-Methyl-N-propyl-pyrrolidinium (P13) Cation and Bis(trifluoromethanesulfonyl)amide and Bis(fluorosulfonyl)amide Anions and Their Binary Systems Including Lithium Salts. J. Chem. Phys. 2010, 133 (19), 194505. (33) Kunze, M.; Jeong, S.; Appetecchi, G. B.; Schönhoff, M.; Winter, M.; Passerini, S. Mixtures of Ionic Liquids for Low Temperature Electrolytes. Electrochim. Acta 2012, 82 (0), 69−74. (34) Kunze, M.; Montanino, M.; Appetecchi, G. B.; Jeong, S.; Schönhoff, M.; Winter, M.; Passerini, S. Melting Behavior and Ionic Conductivity in Hydrophobic Ionic Liquids. J. Phys. Chem. A 2010, 114 (4), 1776−1782. (35) Matsumoto, H.; Sakaebe, H.; Tatsumi, K.; Kikuta, M.; Ishiko, E.; Kono, M. Fast Cycling of Li/LiCoO2 Cell with Low-Viscosity Ionic Liquids Based on Bis(fluorosulfonyl)imide [FSI]−. J. Power Sources 2006, 160 (2), 1308−1313. (36) Moreno, M.; Montanino, M.; Carewska, M.; Appetecchi, G. B.; Jeremias, S.; Passerini, S. Water-Soluble, Triflate-Based, Pyrrolidinium Ionic Liquids. Electrochim. Acta 2013, 99 (0), 108−116. (37) Simons, T. J.; Bayley, P. M.; Zhang, Z.; Howlett, P. C.; MacFarlane, D. R.; Madsen, L. A.; Forsyth, M. Influence of Zn2+ and Water on the Transport Properties of a Pyrrolidinium Dicyanamide Ionic Liquid. J. Phys. Chem. B 2014, 118 (18), 4895−4905. (38) Tokuda, H.; Ishii, K.; Susan, M. A. B. H.; Tsuzuki, S.; Hayamizu, K.; Watanabe, M. Physicochemical Properties and Structures of RoomTemperature Ionic Liquids. 3. Variation of Cationic Structures. J. Phys. Chem. B 2006, 110 (6), 2833−2839. (39) Wang, Y.; Zaghib, K.; Guerfi, A.; Bazito, F. F. C.; Torresi, R. M.; Dahn, J. R. Accelerating Rate Calorimetry Studies of the Reactions between Ionic Liquids and Charged Lithium Ion Battery Electrode Materials. Electrochim. Acta 2007, 52 (22), 6346−6352. (40) Yamaguchi, T.; Mikawa, K.-i.; Koda, S.; Serizawa, N.; Seki, S.; Fujii, K.; Umebayashi, Y. Effects of Lithium Salts on Shear Relaxation 6249

DOI: 10.1021/jp5122678 J. Phys. Chem. B 2015, 119, 6242−6249