J. Phys. Chem. B 2008, 112, 8129–8133
8129
Shear Viscosity of the Ionic Liquid 1-n-Butyl 3-Methylimidazolium Hexafluorophosphate [bmim][PF6] Computed by Reverse Nonequilibrium Molecular Dynamics Wei Zhao,*,† Fre´de´ric Leroy,† Sundaram Balasubramanian,‡ and Florian Mu¨ller-Plathe† Physikalische Chemie, Technische UniVersita¨t Darmstadt, Petersenstrasse 20, D-64287, Germany, and Chemistry and Physics of Materials Unit, Jawaharlal Nehru Centre for AdVanced Scientific Research, Jakkur, Bangalore 560 064, India ReceiVed: February 29, 2008; ReVised Manuscript ReceiVed: April 9, 2008
Reverse nonequilibrium molecular dynamics and equilibrium molecular dynamics simulations were carried out to compute the shear viscosity of the pure ionic liquid system [bmim][PF6] at 300 K. The two methods yielded consistent results which were also compared to experiments. The results showed that the reverse nonequilibrium molecular dynamics (RNEMD) methodology can successfully be applied to computation of highly viscous ionic liquids. Moreover, this study provides a validation of the atomistic force-field developed by Bhargava and Balasubramanian ( J. Chem. Phys. 2007, 127, 114510) for dynamic properties. Introduction The development of ionic liquids has attracted considerable scientific and industrial attention during the past decade.1–3 Because of their unique properties such as low melting point, near-zero volatility, high thermal and electrochemical stability, and electrical conductivity, ionic liquids are considered to be ideal replacements of industrially used organic solvents, which are often not environmentally benign and which are less efficient in production as well as in recycling of materials. To date, a great number of scientific and industrial applications of different ionic liquids have been realized in various fields with an example being the BASIL process whose use has led to an 80 000-fold increase in productivity of the alkoxyphenylphosphine formation compared with the conventional process.1 There are, however, still impediments to further industrial use of ionic liquids, and one impediment is their high viscosity. On the one hand, the high viscosity of ionic liquids leads to difficulties in handling them practically as in processes of pumping, mixing, or stirring. Furthermore, since momentum transport is coupled to other transport properties, a high viscosity goes along with a decrease of mass transfer, heat transfer, and electrical conductivity. Thus, applications of ionic liquids as electrolyte materials in lithium-ion batteries, for instance, are precluded up to now. On the other hand, a low viscosity may cause emulsion formation.4 It is therefore of industrial relevance to design ionic liquids whose viscosity is under control. It has been suggested that the viscosity of ionic liquids is determined by their molecular properties such as their tendency to form hydrogen bonds.2 However, a detailed understanding of the relation between structural and dynamic properties of ionic liquids is still missing, hampering the rational design of low-viscosity products. In recent years, molecular dynamics (MD) simulations5–25 have been performed to obtain fundamental understanding of various properties of ionic liquids in * To whom correspondence should be addressed. Phone: +49-6151-166537. Fax: +49-6151-16-6526. E-mail:
[email protected]. † Technische Universita ¨ t Darmstadt. ‡ Jawaharlal Nehru Centre for Advanced Scientific Research.
the liquid bulk phase or at interfaces. MD methods have also been employed to compute the shear viscosity of different ionic liquids.26–32 The viscosity can be computed by equilibrium MD simulations using the Green-Kubo correlation function formalism or the equivalent Einstein relations30,33 or by analysis of momentum fluctuations.34 However, as Ko¨dderman et al.32 recently suggested, equilibrium MD may be inefficient when applied to highly viscous systems. As an alternative, the system can be driven out of equilibrium by nonequilibrium molecular dynamics methods, and its response can be recorded. The SLLOD algorithm of Evans35 was devised to achieve such a task by measuring the momentum flux (stress) in response to a flow velocity gradient (shear field). However, its application to charged systems might lead to delicate convergence.31 The socalled reverse nonequilibrium molecular dynamics (RNEMD) algorithms36 follow an alternative philosophy by imposing fluxes and measuring the resulting gradients in the system. Thus, in the case of viscosity, a momentum flux is imposed and a velocity gradient is measured. This method has already been employed to compute the shear viscosity of several liquid systems including both simple and polymeric liquids.36–39 Recently, Kelkar and Maginn31 successfully applied it to the aqueous [emim][Tf2N] ionic liquid system over a range of temperatures and compositions. Still, the largest viscosity computed in that work was about 40 cP. In the present work, we demonstrate the ability of the RNEMD method to study bulk liquid 1-nbutyl 3-methylimidazolium hexafluorophosphate [bmim][PF6](see Figure 1 for chemical structure) which is of much higher viscosity. The present paper is organized as follows: in the first part, some details about the atomistic force-field and the simulations are given; in the second part, the RNEMD simulations data are presented and discussed; finally, some conclusions are drawn. Model Description and Methods All atoms of the system are represented in the force field of [bmim][PF6] used here. This model also takes into account all the molecular internal degrees of freedom (bonds, bond angles, and torsional angles) as well as nonbonded repulsion-dispersion
10.1021/jp8017869 CCC: $40.75 2008 American Chemical Society Published on Web 06/18/2008
8130 J. Phys. Chem. B, Vol. 112, No. 27, 2008
Zhao et al.
Figure 1. Chemical structures of 1-n-butyl-3-methylimidazolium ([bmim]; bottom) and hexafluorophosphate (PF6).
and electrostatic interactions such that the total interaction potential writes bonds
Utotal )
∑ ij
dihedrals
∑ ijkl
kr,ij (r - r0,ij)2 + 2 ij
angles
∑ ijk
kφ,ijk (φijk - φ0,ijk)2 + 2
kτ,ijkl [1 - cos pijkl(τijkl - τ0,ijkl)] + 2
( )) σij rij
6
EL
+
∑ ij
(
LJ
(( )
∑ 4ε ij
σij rij
cutoff
-
)
(1)
For the description of the terms, see ref 23. The parameters were taken from the model of Bhargava and Balasubramanian.40 A cutoff of 1.2 nm was applied for the nonbonded interactions. The mixed Lennard-Jones parameters were derived from selfparameters using the Lorenz-Berthelot mixing rule.41 Also within a same molecule, 1-2 and 1-3 nonbonded interactions were excluded, and 1-4 nonbonded interaction strengths were divided by a factor of 2 as it is usual in all-atom OPLS-type force fields.42 The reaction-field methodology was employed to compute electrostatic interactions with a dielectric constant εRF of 11.4.43 We here briefly describe the basic ideas of the RNEMD method. More details can be found elsewhere.38 The shear-ratedependent shear viscosity (γ˙ ) of fluids is defined as the proportionality constant between a shear field γ˙ ) ∂Vx/∂z and a flux of transverse linear momentum jz(px) when linear response holds
jz(px) ) -η(γ˘ )
∂Vx ∂z
momentum flux depends on the total momentum ∆p transferred in the z direction per time ∆t and cross-sectional area Lx × Ly
jz(px) )
12
εRF - 1 rij2 qiqj 1 + 4πε0 rij 2εRF + 1 r3
Figure 2. Temperature profiles in systems at different swap periods. The temperature in each slab is calculated using only the random part of the total kinetic energy of the slab by subtracting the steady flow in the x direction. Symmetry-equivalent data are averaged; only one-half of the simulation box is shown.
(2)
The flux jz(px) is the x component of momentum px that is transported in the z direction per unit time and unit area. The direction chosen for flux and momentum projections can be chosen arbitrary as long as they are perpendicular. The RNEMD method imposes an unphysical momentum transfer within the simulation box in the z direction. To this end, the simulation box is divided into N virtual slabs. Periodically, the center-of-mass velocities of two molecules are exchanged: the first one has the largest negative x component of the center-of-mass momentum in slab 1 and the other one belonging to the same molecular species has the largest positive x component of the center-of-mass momentum in slab N/2 + 1. As a response, a physical momentum flow arises in the z direction, and a linear velocity profile appears. The imposed
∆p 2LxLy∆t
(3)
After a transition period, the system converges toward a steady state characterized by a constant gradient in transverse velocity (∂Vx/∂z). The later part of simulation trajectory is used for accumulation of the velocity gradient and the momentum transfer. By computing the slope of the velocity profile and using eq 2, one obtains the shear viscosity at a given shear rate. The RNEMD simulations were performed at 300 K with different shear rates. The simulation boxes contained 768 ion pairs (24 576 atoms) and were divided into 20 slabs. The initial configurations were prepared from pre-equilibrated cubic simulation boxes containing 256 ion pairs by replicating them three times in the z-direction, which resulted in L × L × 3L simulation boxes (L ) 4.455 nm). At 300 K, the velocity exchange periods were taken to be 5, 10, 20, 40, 80, 100, 120, and 160 time steps. The simulations of systems at swap periods of 5, 10, and 20 were performed for 14 ns, and the last 8 ns were used for data analysis. For systems at swap periods of 40, 60, 80, 100, 120, and 160, time steps were simulated for at least 32 ns and at least 15 ns were used for data analysis. In all simulations, the temperature was kept constant using the Berendsen thermostat44 with a coupling constant of 1.0 ps. All simulations were performed using our molecular simulation package YASP.45,46 To provide a comparison with our results obtained in RNEMD simulations, we also carried out an equilibrium molecular dynamics simulation of a [bmim][PF6] system containing 256 ion pairs for 90 ns of which the last 80 ns were used for data analysis. Each simulation was run in parallel on an IBM p375 cluster. In total, the simulations took ∼60 000 h of CPU time. Results and Discussions 1. Temperature and Density Profiles. In reverse nonequilibrium molecular dynamics, undirected motion (heat) is turned into directed motion (shear flow) by exchanging of momentum of selected molecules in two regimes. This leads to the phenomenon that the temperature is not uniform throughout the simulation box. Figure 2 shows the temperature profiles of systems at different exchange intervals. At strong perturbations for swap period W ) 5, 10, 20, and 40 time steps, the temperature profiles exhibit a parabola-like feature with differ-
Viscosity of [bmim][PF6] by Molecular Dynamics
Figure 3. Density profiles of systems at swap periods of 5, 10, 20, and 80. Dash-dotted line denotes the equilibrated density in absence of shear. Symmetry-equivalent data are averaged; only one-half of the simulation box is shown.
ences in temperature of 63.1, 40.7, 10.6, and 1.9 K between hot and cold regimes. This parabola-like behavior has been found for other molecular fluids in our previous work.36 As the swap period goes beyond 40 time steps, the temperature profiles become less well-defined, as shown in Figure 2b, because of the dominance of statistical fluctuations. The temperature differences within the simulation box are very small. For example, at an exchange period of 100, the temperature variation is less than 1.3 K, which is within the statistical fluctuations of about 7 K. Figure 3 shows the density profiles in these systems at different shear rates. It is not surprising that the density profiles are inversely coupled to the temperature profiles. At exchange period of 5, for example, the density in the cold region (e.g., exchange slabs) is 5.5% higher than in the hot areas. As the exchange period becomes 20, this variation drops to around 2%, which is within the density fluctuations of 3% in these regions. For weak perturbations (exchange period >40), the residual small temperature variation does not cause a systematic and significant density variation. On the basis of the observations about the temperature and density profiles, we can safely assume that the weaker perturbations (swap period of 40 time steps or larger) leave the system in a uniform phase. 2. Velocity Profiles. In the present simulations, velocity profiles were obtained by calculating the velocity in the x direction averaged over all ions in one slab and over the trajectory. The gradient of the velocity profile represents the strength of the shear field. A good linearity of the velocity profile is required to calculate the shear viscosity, and one has to choose carefully the shear rates at which linear response holds. In all simulations, [bmim] cations were taken as the species for exchanging momentum on the basis of an observation that the choice of [PF6] anion might lead to severe nonlinearity in the velocity profiles of which the mechanism needs to be understood through further investigation. Our analysis also shows that no electric current was induced by exchanging the momentum of cations only. For the present ionic liquid, nonlinear velocity profiles are found for strong shear rates (swap periods of 5 and 10 time steps, Figure 4). As the swap periods go beyond 20 time steps, the velocity profiles turn out to be linear and can be used to calculate the shear viscosity using eq 2. At weaker perturbations (Figure 5), the velocity gradients are linear but take longer to converge. However, perturbations that are too weak (at an exchange period of 200 time steps, for example) are found to exhibit low signal-to-noise ratios even after a long computing time. 3. Shear Viscosity. For Newtonian molecular fluids of low viscosity, the viscosity is independent of the shear rate for a
J. Phys. Chem. B, Vol. 112, No. 27, 2008 8131
Figure 4. Velocity profiles in the x direction for exchange frequencies 5, 10, 20, and 40. Symmetry-equivalent data are averaged; only onehalf of the simulation box is shown.
Figure 5. Velocity profiles in the x direction for exchange frequencies 80 (solid), 100 (long dashed), 120 (short dashed), and 160 (dotted). The straight lines indicate least-squares fits to the data sets. Symmetryequivalent data are averaged; only one-half of the simulation box is shown.
large range of shear rates. Only for very high shear rates do nonlinearities such as shear-thinning appear. One can, therefore, choose different shear rates within the plateau range. The present system, however, does not have such a plateau, and a dependence of the shear viscosity on the chosen shear rate is found. The shear viscosity, at a swap period of 20 time steps, for example, is 12.4 ( 0.2 cP, which is 1 order of magnitude lower than the value of 118.8 ( 31.1 cP at a swap period of 160. The fact that the shear viscosity is smaller at stronger perturbation (shear thinning) probably results from the structural alterations under strong shear flow, such as reorientation of [bmim] cations and fewer hydrogen bonds. A more detailed study will be reported later. Extrapolating the zero-shear viscosity from a nonconstant shear viscosity requires a model. Although mode-coupling theory has been used as an extrapolation procedure by several authors,31,47 these authors have also noticed that it can only be applied to a restricted shear-rate domain. Alternatively, one can use the empirical three-parameter Carreau equation31 to fit the viscosity as a function of the shear rate
η)
η0 [1 + (λγ˘ )2]p
(4)
where η0 is the Newtonian (zero-shear) viscosity, λ is a characteristic time constant, and p is an exponent. λ approximately corresponds to the rotational time of the cation.31 Using this procedure (see Figure 6), we obtained a value of 139.6 ( 9.7 cP (with λ of 9.0 ns and p of 0.28) as an estimate of the zero-shear viscosity of the present system. This value is lower than the 185.2 cP48 and 228.8 cP49 at 300 K observed experimentally (data obtained from VFT fitting of experimental observations). This underestimation is in line with faster ion dynamics described in the present model. The diffusion coef-
8132 J. Phys. Chem. B, Vol. 112, No. 27, 2008
Zhao et al. Conclusions
Figure 6. Shear viscosities at different swap periods of 40, 60, 80, 100, 120, and 160. Solid line represents the fit to the Carreau equation. The calculation of the error bar of the shear viscosity has been described in ref 38.
Figure 7. Mean square displacements of [bmim] cations, [PF6] anions, and all ions on average.
ficients of cations and anions were calculated to be 1.29 × 10-11 m2 s-1and 9.0 × 10-12 m2 s-1 (see Figure 7 for the mean square displacements of cations and anions), which are about 60% higher than the experimental values of 8.0 × 10-12 m2 s-1 and 5.9 × 10-12 m2 s-1, respectively.49 An estimate of the zero-shear viscosity was also computed by equilibrium MD using the Einstein expression. It has the practical advantage that it does not require the frequent output of the stress tensor components to take into account the fast oscillations in the Green-Kubo integral.33 The zero-shear viscosity is computed as
η)
d V lim 20kBT tf∞ dt
[∑ 〈 R)x,y,z
[PRR(t) - PRR(0)]2 〉 + 2
∑ 〈[PRβ(t) - PRβ(0)]2〉] (5)
β>R
which appears as the asymptotic time derivative of the generalized mean-squared displacement of the following quantity:
PRβ(t) )
∫0t ΠRβ(s)ds
(6)
where
1 1 ΠRβ(t) ) [σRβ(t) + σβR(t)] - δRβ trace[σ(t)] 2 3
(7)
σ(t) is the instantaneous stress tensor including the kinetic contribution Ekin and the virial contribution Ξ
2 σ ) (E - Ξ) V
(8)
A trajectory of 80 ns was used, and the time derivative appearing in eq 5 was computed with a 40 ns correlation length. The zeroshear viscosity computed in this way converges after about 10.5 ns and is found to be 127 ( 17 cP, which in good agreement with the RNEMD result.
In the present work, the reverse nonequilibrium molecular dynamics approach has been demonstrated to be able to study the shear viscosity of an ionic liquid of high viscosity. For the studied system [bmim][PF6], our results are in excellent agreement with computationally more demanding equilibrium MD simulations. The force field used produces a good agreement with experiment, although the calculated viscosity appears to be below the experimental value by 25-40% depending on which experiment is used as reference. This is, however, a vast improvement over some previous force fields which typically overestimated the viscosities of ionic liquids by more than 1 order of magnitude. In our study, it was found that the choice of [PF6] anions for exchanging momentum might lead to severe nonlinearity in velocity profile. The mechanism of this unexpected behavior is still not clear and will be investigated in our future work. Also, further studies are underway to elucidate the effects of shear on the structural and dynamic properties of [bmim]-based ionic liquid systems. Acknowledgment. This work is supported by the Deutsche Forschungsgemeinschaft (DFG) by the SPP1191 “Ionic liquids” priority program. We thank Deutscher Akademischer Austauschdienst (DAAD) for providing travel grants. We thank HHLR supercomputer centre at the Technische Universita¨t Darmstadt for providing computing resources. S.B. thanks Department of Science and Technology of India (DST) for support. We also gratefully thank Prof. P. Wasserscheid (Erlangen, Germany), Prof. A. Heintz (Rostock, Germany), Prof. B. Kirchner (Leipzig, Germany), and Prof. H. Weinga¨rtner (Bochum, Germany) for fruitful discussions. References and Notes (1) Rogers, R. D.; Seddon, K. R. Science 2004, 302, 792–793. (2) Wasserscheid, P.; Keim, W. Angew. Chem., Int. Ed. 2000, 39, 3772– 3789. (3) Weinga¨rtner, H. Angew. Chem., Int. Ed. 2008, 47, 654–670. (4) Anjan, S. T. Chem. Eng. Prog. 2006, 102, 30–39. (5) Zhao, W.; Eslami, H.; Cavalcanti, W. L.; Mu¨ller-Plathe, F. Z. Phys. Chem. 2007, 221, 1647–1662. (6) Cadena, C.; Anthony, J. L.; Shah, J. K.; Morrow, T. I.; Brennecke, J. F.; Maginn, E. J. J. Am. Chem. Soc. 2003, 126, 5300–5308. (7) Cadena, C.; Maginn, E. J. J. Phys. Chem. B 2006, 110, 18026– 18039. (8) Cadena, C.; Zhao, Q.; Snurr, R. Q.; Maginn, E. J. J. Phys. Chem. B 2006, 110, 2821. (9) Canongia Lopes, J. N.; Gomes, M. F. C.; Pa´dua, A. A. H. J. Phys. Chem. B 2006, 110, 16816–16818. (10) Canongia Lopes, J. N.; Pa´dua, A. A. H. J. Phys. Chem. B 2006, 110, 19586–19592. (11) Canongia Lopes, J. N.; Pa´dua, A. A. H. J. Phys. Chem. B 2004, 108, 16893–16898. (12) Canongia Lopes, J. N.; Deschamps, J.; Pa´dua, A. A. H. J. Phys. Chem. B 2004, 108, 11250. (13) Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. A. H. J. Phys. Chem. B 2004, 108, 2038–2047. (14) Pa´dua, A. A. H.; Gomes, M. F. C.; Canongia Lopes, J. N. Acc. Chem. Res. 2007, 40, 1087–1096. (15) Lynden-Bell, R. M.; Youngs, T. G. A. Mol. Simul. 2006, 32, 1025– 1033. (16) Lynden-Bell, R. M.; del Popolo, M. G.; Youngs, T. G. A.; Kohanoff, J.; Hanke, C. G. Acc. Chem. Res. 2007, 40, 1138–1145. (17) Lynden-Bell, R. M.; del Popolo, M. G. Phys. Chem. Chem. Phys. 2006, 8, 949–954. (18) Lynden-Bell, R. M.; Kohanoff, J.; del Popolo, M. G. Faraday Discuss. 2005, 129, 57–67. (19) Lynden-Bell, R. M.; del Popolo, M. G. Phys. Chem. Chem. Phys. 2005, 8, 949–954. (20) Hanke, C. G.; Lynden-Bell, R. M. J. Phys. Chem. B 2003, 107, 10873–10878.
Viscosity of [bmim][PF6] by Molecular Dynamics (21) Bhargava, B. L.; Klein, M. L.; Balasubramanian, S. Chem. Phys. Chem. 2008, 9, 67. (22) Bhargava, B. L.; Balasubramanian, S. J. Am. Chem. Soc. 2006, 128, 10073–10078. (23) Sloutskin, E.; Lynden-Bell, R. M.; Balasubramanian, S.; Deutsch, M. J. Chem. Phys. 2006, 125, 174715. (24) Hanke, C. G.; Price, S. L.; Lynden-Bell, R. M. Mol. Phys. 2001, 99, 801–809. (25) Huddleston, J. G.; Visser, A. E.; Reichert, W. M.; Willauer, H. D.; Broker, G. A.; Rogers, R. D. Green Chem. 2001, 3, 156–164. (26) Yan, T.; J., B. C.; del Popolo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 11877–11881. (27) Bhargava, B. L.; Balasubramanian, S. J. Chem. Phys. 2005, 123, 144505. (28) Micaelo, N. M.; Baptista, A. M.; M., S. C. J. Phys. Chem. B 2006, 110, 14444–14451. (29) Rey-Castro, C.; Vega, L. F. J. Phys. Chem. B 2006, 110, 14426– 14435. (30) Borodin, O.; Smith, G. D. J. Phys. Chem. B 2006, 110, 11481– 11490. (31) Kelkar, M. S.; Maginn, E. J. J. Phys. Chem. B 2007, 111, 4867– 4876. (32) Ko¨ddermann, T.; Paschek, D.; Ludwig, R. ChemPhysChem. 2007, 8, 2264–2470. (33) Mondello, M.; Grest, G. S. J. Chem. Phys. 1997, 106, 9327–9336. (34) Hess, B. J. Chem. Phys. 2002, 116, 209–217.
J. Phys. Chem. B, Vol. 112, No. 27, 2008 8133 (35) Evans, D. J. Phys. ReV. A 1991, 44, 3630–3632. (36) Mu¨ller-Plathe, F. Phys. ReV. E 1999, 59, 4894–4899. (37) Cavalcanti, W. L.; Chen, X.; Mu¨ller-Plathe, F. Phys. Status Solidi 2007, 204, 935–939. (38) Bordat, P.; Mu¨ller-Plathe, F. J. Chem. Phys. 2002, 116, 3362–3369. (39) Chen, X.; Carbone, P.; Cavalcanti, W. L.; Milano, P.; Mu¨ller-Plathe, F. Macromolecules 2007, 40, 8087–8095. (40) Bhargava, B. L.; Balasubramanian, S. J. Chem. Phys. 2007, 127, 114510. (41) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press: Oxford, U.K., 1989. (42) Jogensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236. (43) Weinga¨rtner, H. Z. Phys. Chem 2006, 220, 1395–1405. (44) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (45) Mu¨ller-Plathe, F. Comput. Phys. Commun. 1993, 78, 77–94. (46) Tarmyshov, K.; Mu¨ller-Plathe, F. J. Chem. Inf. Mod. 2005, 45, 1943. (47) Travis, K.; Searles, D. J.; Evans, D. J. J. Mol. Phys. 1998, 95, 195. (48) Baker, S. N.; Baker, G. A.; Kane, M. A.; Bright, F. V. J. Phys. Chem. B 2001, 105, 9663. (49) Tokuda, H.; Hayamizu, K.; Ishii, K.; Susan, A. B. H.; Watanabe, M. J. Phys. Chem. B 2004, 108, 16593.
JP8017869
