8954
J. Phys. Chem. B 2010, 114, 8954–8960
Calculating Thermodynamic Properties of an Ionic Liquid with Monte Carlo Simulations with an Orthorhombic and a Cubic Simulation Box Bjo¨rn Wittich* and Ulrich K. Deiters Department of Physical Chemistry, UniVersity of Cologne, Cologne, Germany ReceiVed: February 24, 2010; ReVised Manuscript ReceiVed: May 28, 2010
In contrast to the common usage of cubic simulation boxes, in this work simulations of the ionic liquid 1-n-butyl-3-methyl-imidazolium hexafluorophosphate ([bmim][PF6]) were carried out in a dynamic orthorhombic simulation box over a temperature range from 313 to 373 K in a canonical harmonical simulation ensemble (NpT) with a united-atom potential based on quantum chemistry. The solubilities of the gases CO2, CO, H2, O2, C2H4, and H2O at infinite dilution were determined by means of the Widom test particle method; the results are compared with experimental data and simulation results obtained with a cubic simulation box. For gas potentials containing partial charges the results are in good agreement with the experimental data. Introduction Ionic liquids (ILs) are salts that are liquid at temperatures below 100 °C. In the previous decades, these compounds have gained importance because of their unusual properties. Ionic liquids are excellent solvents for polar compounds because of their polar nature. They have a very low vapor pressure, so that separation of solutes by distillation is relatively easy. Other important properties are the enhanced oxidation, reduction, and thermal stability. Some years ago the first industrial processes utilizing ILs were announced. An important thermophysical property is the solubility of gases in ionic liquids. However, it is often difficult to determine this solubility experimentally, and sometimes the error margins of the experimental data are of the order of magnitude of the values themselves. In such cases, computer simulations can be an alternative way to confirm and to predict gas solubilities. One of the first simulations of gas solubilities in ionic liquids was carried out by Shah and Maginn.1 They developed a unitedatom model of [bmim][PF6] based on quantum chemistry and performed NpT Monte Carlo simulations. With the obtained equilibrium configurations the residual chemical potential of carbon dioxide was calculated by means of the Widom test particle method. It turned out that the calculated Henry’s constants were too small by a factor of 3, indicating a too large solubility. A different approach was used by Urukova et al.2 Here Gibbs ensemble simulations were performed to calculate the solubilities of CO2, CO, and H2 in [bmim][PF6] and then extrapolated to zero solute concentration. The first presented results of Urukova et al. stood in excellent agreement with experimental data,2 but a few months later the results were revoked because of an conversion mistake.3 After correcting this mistake, the results turned out to be nearly the same as those of Shah and Maginn. In a previous paper,4 we used lithium iodide as an example to demonstrate how the geometry of the simulation box may influence the outcome of a computer simulation in solid or dense fluid phases. Now the application temperatures of ionic liquids are usually relatively close to the melting point, the density is high, and so the restriction of the simulation box to a cubic shape might influence the outcome of simulations for common system sizes of 250 to 500 particles. In this paper, we compare simulations of [bmim][PF6] using orthorhombic boxes to * To whom correspondence should be addressed: E-mail: Bjoern.Wittich@ uni-koeln.de. Phone+49 (0)221 470-4543. Fax+49 (0)221 470-4900.
simulations with cubic boxes. Structural information, internal energy, and Henry’s constants are reported. Simulation Details NpT Simulation Methodology. An isothermal-isobaric canonical harmonical ensemble was employed to obtain the equilibrium configurations of the pure ionic liquid. The start configuration was generated by randomly placing 125 cations and 125 anions into a huge simulation box. The preequilibration was carried out in the cubic simulation box to prevent the simulation box distorting too much in this early stage. When the volume of the simulation box had converged, a system configuration was selected and used for two different simulations. One was performed with a cubic simulation box again, in the other simulation the dimensions of the simulation box could change independently, so that an orthorhombic simulation box resulted. In simulations of solid phases, we showed that the difference between orthorhombic and parallelepiped cells is negligible, hence no attempt was made here to let the box angles vary.4 In simulations of isotropic phases, as in this work, the orthorhombic box never assumed a thin and elongated shape, but remained roughly cubic as shown in Figure 1. Test simulations with different particle numbers were performed in order to ascertain that the corrections for the truncation of shortand long-ranged potentials were adequate. Equilibration usually took about 105 cycles; one cycle consists of N particle movements and three attempts to change the length of a randomly chosen box side. The particle movements consisted of linear displacements, followed by rotation about their centers of gravity. The shifting and rotating parameters were adjusted for each species to yield an acceptance ratio of 0.3. The acceptance criterion for the particle movement (pm) is
ppm(o f n) )
{
if ∆Uofn < 0 1, exp(-β(∆Uofn)), otherwise
(1)
The criterion for the volume change (vc) is also based on the Metropolis acceptance scheme
pvc(o f n) )
{
1, if a > 1 a, otherwise
(2)
where a ) exp (-β[∆U + p∆V - (N + 1)β-1 ln(Vn/Vo)]).
10.1021/jp101676m 2010 American Chemical Society Published on Web 06/21/2010
Properties of an Ionic Liquid with Monte Carlo Simulations
J. Phys. Chem. B, Vol. 114, No. 27, 2010 8955
Figure 1. Fluctuation of the orthorhombic simulation box in the equilibrium. black line, x-dimension; red line, y-dimension; blue line, z-dimension.
Figure 2. Illustration of the ions (a) [bmim]+ cation (b) [PF6]- anion used in this work (see Table 1).
The configuration obtained in a previous simulation run was taken as start configuration for the new simulation run at a slightly different temperature. We found however that it is very inefficient to use a configuration of a simulation run at lower temperature as initial configuration of a new run at a higher temperature, because equilibration takes a very long time. At low temperatures, the disorder of particles and also the amount of vacancies in the system is reduced. It takes a very long time to create these vacancies again when the system is heated; the system shows a kind of hysteresis. Molecular Models. The pair potential of the ionic liquid [bmim][PF6] is an approximate united-atom model (UA2) based on quantum chemistry proposed by Shah et al.1 In this model all hydrogen atoms of the cation are incorporated into the carbon atoms to which they are bound. This reduces the number of interaction sites from 26 to 11. Urukova et al. used a similar model2 of Shah and Maginn (UA1) where also the fluorine atoms are incorporated in the phosphorus central atom. Although Shah and Maginn considered a term for intramolecular rotations, these contributions where omitted in this work in accordance with the work of Urukova et al. where these contributions were also neglected. The potential models and their parameters for the dissolved gases are given in Figures 2-9. For the production runs of carbon dioxide the EPM2-potential of Harris and Yung5 and
Figure 3. Potential of CO2.5
for source code testing the TraPPe-potential of Potoff and Siepmann6 were used. Both results are in good agreement. Hydrogen was described with the model of Cracknell,7 which consists of two Lennard-Jones centers. A similar model was used for oxygen, where a potential of Miyano8 was used. For carbon monoxide, two different models were used, namely the simplified model of Bohn et al.,9 which is a 2-center LJ potential, and the more accurate model of Straub and Karplus10 where the electrical partial charges are considered explicitly. The water model of Berendsen et al.11 describes water by three partial charges and one LJ center (located at the oxygen). Finally ethene is represented by two LJ centers and three partial charges. The charges are located at the C atoms and the center of gravity.
8956
J. Phys. Chem. B, Vol. 114, No. 27, 2010
Wittich and Deiters TABLE 1: Cartesian Coordinates of the UA Model of [bmim][PF6]1a
Figure 4. Potential of CO (with ionic charges).10
Figure 5. Potential of CO (LJ model).9
atom
X
Y
Z
σii/Å
εii/kJ/mol
q/10-19C
N1 N3 C5 C7 C9 P10 F2 F4 F6 C2 C4 C6 C8 C10 F1 F3 F5
-2.048 -2.041 -3.279 -1.574 -1.340 1.995 0.502 2.462 1.366 -1.322 -3.274 -1.542 -1.993 0.173 1.428 2.529 3.400
0.682 -1.177 0.082 1.940 3.348 -1.001 -1.550 -0.861 0.553 -0.100 -1.083 -2.310 3.187 3.582 -1.098 -2.527 -0.413
0.272 -0.852 0.469 0.888 -1.278 0.512 1.005 2.055 0.557 -0.530 -0.232 -1.640 0.103 -1.238 -1.066 0.414 -0.031
3.250 3.250 3.880 3.905 3.905 3.740 3.118 3.118 3.118 3.880 3.880 3.775 3.905 3.905 3.118 3.118 3.118
0.710 0.710 0.443 0.493 0.493 0.836 0.255 0.255 0.255 0.443 0.443 0.865 0.493 0.732 0.255 0.255 0.255
0.111 0.133 -0.010 0.195 0.128 1.460 -0.394 -0.394 -0.394 0.233 0.040 0.183 -0.066 -0.043 -0.394 -0.394 -0.394
a
8
Figure 6. Potential of O2.
Figure 7. Potential of H2.7
The associated molecules are shown in Figure 2.
particle numbers, we found that the truncation error is very small and no correction is necessary. For the Lennard-Jones interactions between unlike sites, standard Lorentz-Berthelot combining rules were applied. The cross size parameter is determined by the arithmetic mean and the cross energy well depth is determined by the geometric mean of the parameters of the pure components. Calculating Henry’s Constant with the Widom Test Particle Method. To obtain Henry’s constant from molecular simulations, it is necessary to determine the residual chemical potential of the dissolved gas. The solubility of a small amount of gas (s, solute) in the IL (l, liquid) can be described by Henry’s constant13
( )
kH ) lim FlkBT exp xsf0
Figure 8. Potential of H2O.11
µsres,l kBT
(3)
The residual chemical potential can be estimated with the Widom test particle method, which has been described elsewhere.14,15 It should be mentioned that the Ewald sum algorithm had to be modified because the inserted test particle does not occur in the periodic images as assumed in the original Ewald formulation. Simulation Results
Figure 9. Potential of C2H4.12
Interaction Potentials. The minimum image convention was applied to the calculation of the short ranged Lennard-Jones (LJ) interactions, whereas the electrostatic contribution of the partial charges was treated with the Ewald sum method, after implementing the necessary modifications for noncubic simulation boxes. For the short-range interactions, the cutoff distance is g9 Å. Because of the long-range correlations of this ionic system, it is not possible to compute a tail correction. But by computing the short-range energies of systems with different
The simulation ensemble consisted of 125 cation-anion pairs at temperatures 315.15 or 373.15 K. At both temperatures, simulations in a cubic and an orthorhombic simulation box were carried out at a pressure of 1 bar. The usual length of a simulation was 105 cycles. The density data and the data of Henry’s constant of carbon dioxide at 343.15 K were obtained with an ensemble of 64 cation-anion pairs. With this smaller ensemble, we also simulated temperatures of 313.15 and 373.15 K, but the resulting values of Henry’s constant were exactly the same as in the larger 125 cation-anion-pairs ensemble, so there is no need to present these data here. Density. The density as an average of the single system configurations is given in Table 2 as well as Figure 10 and compared with experimental data. These data agree well with the simulation data of Shah and Maginn,1 who also found a lower density in comparison to experimental data with this molecular potential.
Properties of an Ionic Liquid with Monte Carlo Simulations TABLE 2: Density of [bmim][PF6] F/g cm-3 T/K
simul orthorh
simul cubic
exp.16
313.15 343.15 373.15
1.255 ( 0.008 1.270 ( 0.010 1.293 ( 0.007
1.260 ( 0.010
1.3489 ( 0.0080 1.3239 ( 0.0078
1.290 ( 0.007
Radial Distribution Function. The radial cation-anion distribution function (RDF) of the center of mass is shown in Figure 11. The positions of the maxima and minima are very similar to those observed by Shah and Maginn1 but the g(r) value for the RDF of mass centers is larger than observed there. Orientation Autocorrelation Function. To analyze whether the molecules are able to rotate during the simulation run, the orientation autocorrelation function was determined for each species. The orientation of the cation is described by a unit vector b a that is perpendicular to the imidazole ring, and the
J. Phys. Chem. B, Vol. 114, No. 27, 2010 8957 orientation of the anion is specified by a unit vector pointing through two opposite fluorine atoms of the hexafluorophosphate anion. In Figure 12, the orientation autocorrelation function is shown as a function of the number of simulation cycles. Its definition is given by eq 4
cAA(t) )
CAA(t) a 0〉 〈a bt · b ) CAA(0) 〈a b0 · b a 0〉
(4)
It should be noticed that the orientation correlation of the anions vanished after a few 100 cycles (see Figure 13). This can be explained by the almost globular shape of the anion. The rotation barrier in the temperature range of this work is so small that a rotation of the anion is smoothly possible. The behavior of the cation was quite different. Even after a very long simulation time, the autocorrelation function decreased to
Figure 10. Density plot versus temperature. Black line, exp. data; red line, simulated data.
Figure 11. The center-of-mass radial distribution function between cation and anion at 313 K (red line) and 373 K (black line).
8958
J. Phys. Chem. B, Vol. 114, No. 27, 2010
Wittich and Deiters
Figure 12. Orientation autocorrelation function of the cation plotted versus the number of simulation cycles: (black line) 373 K, (blue line) 343 K, (red line) at 313 K. (a) Orthorhombic simulation box and (b) cubic box.
Figure 13. Orientation autocorrelation function of the anion plotted versus the number of simulation cycles: (black line) 373 K, (blue line) 343 K, (red line) at 313 K. (a) Orthorhombic simulation box and (b) cubic simulation box.
about 0.7 only (see Figure 12). That is not really surprising, because the planar geometry of the molecule and its many coordination sites greatly hinder the rotation. Between the cubic and orthorhombic simulation runs there are no significant differences. Orientation Cross Correlation. Since the rotation of the cations is hindered, it is important to check for a correlation between the orientations of the cations. The correlation of the anion orientations is not important, because it is clear that they can rotate freely. At each cycle the cation-cation orientation correlation is characterized with the order parameter S of the Maier-Saupe theory
3 1 S ) 〈cos2 (θ)〉 2 2
(5)
where θ is the angle between the orientation vectors of two cations. A value of 0 indicates that there is no spatial correlation between the molecules, and the phase is isotropic. A value of 1 describes a full alignment of the molecules. The order parameter was calculated at each cycle of simulation, using the previously definied orientation vectors. The results are plotted in Figure 14. All determined values of S are located in an interval between 0.025 and -0.007. This indicates that the distribution is isotropic.
Diffusion. The most interesting quantity is the diffusion of molecules. Although there is no time scale in a Monte Carlo simulation, it is possible to plot the mean square displacement versus the number of simulation cycles. Because of the Einstein-Smulochowski equation, a linear trend of the plot should be expected. As the observed diffusion rates of cations and anions are nearly equal, in Figure 15 the averaged squared particle diffusion is plotted versus the number of cycles. Obviously the mean square displacement shows a dependence on the geometry of the simulation box. In the case of the orthorhombic simulation box the trend is nearly linear. For the cubic simulation box, however, we found a characteristic point of inflection especially at lower temperatures. Futhermore, it should be noted that the “diffusion rate” (squared displacement per number of cycles) is nearly two times larger for the orthorhombic simulation. The differences in the diffusion behavior in the orthorhombic and the cubic simulation box can probably be attributed to the restriction of the configuration room (so-called bottlenecks). It is possible that some energetically favorable configurations could not be left any longer, because all reachable new configurations required too much energy. Figure 15a,b merely show portions of the simulation runs. The whole runs achieved mean-squared displacements of about 6.5 Å, which is comparable to the box length. Comparison with experimental diffusion data17 shows that our simulations spanned
Properties of an Ionic Liquid with Monte Carlo Simulations
J. Phys. Chem. B, Vol. 114, No. 27, 2010 8959
Figure 14. Correlation of the cation-cation orientations at different temperatures: (black line) 373 K, (blue line) 343 K, (red line) at 313 K. (a) Orthorhombic box and (b) cubic box.
Figure 15. Comparison of the averaged mean square displacement at 373 K (black line) and 313 K (red line). (a) Panel shows the case of the orthorhombic box and (b) with a cubic box. The nonlinear plot in the cubic case is clearly demonstrated.
a real time of about 13 ps during the production run. For the equilibration a much longer time span was passed. Henry’s Constant. Henry’s constants of several gases were obtained from the residual chemical potential by means of Widom’s test particle method. Test particles (105) were inserted in 2000 randomly chosen system configurations of the NpT simulation run. This procedure was repeated 5 times for each gas species at every temperature. Henry’s constant was obtained by averaging the single values. Table 3 contains the obtained values. Figure 16 shows the dependency of Henry’s constant on the number of insertion attempts. Energy of the Condensed Phase. During the simulation the difference of the internal energies of the liquid phase and a hypothetical gas phase consisting of ions at infinite distance was obtained from the equilibrated configurations. Although it is not possible to compare to experimental data, the obtained data are presented in Table 4
∆U ) Ugas,inf .dist. - Uliquid
(6)
Discussion MC simulations of the ionic liquid [bmim][PF6], modeled as a united-atom multicenter LJ molecule with partial charges on each site at 313-373 K and 1 bar were performed using a NpT-
TABLE 3: Overview of the Determined Henry’s Constants (mole fraction scale) of Various Gases in [bmim][PF6] gas CO2 CO2 C2 H 4 C2 H 4 COion. COion. COLJ COLJ 21 O2 H2 H2 H 2O H 2O
T/K kH, sim, orthorh./bar 313 373 313 373 313 373 313 373 313 373 313 373 313 373
60 ( 7 150 ( 4 96 ( 18 117 ( 19 2326 ( 166 3284 ( 134 4697 ( 66 5123 ( 46 2654 ( 40 3518 ( 40 6040 ( 15 6351 ( 46 0.07 ( 0.04 1.54 ( 0.16
kH, sim, cub./bar 14 ( 1 109 ( 5 162 ( 4 72 ( 7 422 ( 12 1685 ( 63 1700 ( 28 3062 ( 44 1339 ( 14 2281( 50 4736 ( 27 4790 ( 32 0.09 ( 0.03 1.16 ( 0.54
kH, exp./bar 65.3,18 71.2 ( 2.619 157.418 190 ( 2319 1937 ( 2720 1971 ( 2720 1937 ( 2720 1971 ( 2720 1807 ( 1821,22 1770 ( 4021,22 4154 ( 4523 3059 ( 3923 24 24
ensemble with a variable simulation box shape. The orientation autocorrelation function shows why the ideal crystal structure is inappropriate as a start configuration. The equilibration in the dense phase is hindered because of the high rotational barriers of the planar cation. Our method to prepare the ensemble from an ideal gas state ensured that there were no spatial orientational correlations. The orthorhombic simulation box geometry, developed for noncubic solids, was not expected to have any effect for simulations of isotropic liquids. Indeed, the
8960
J. Phys. Chem. B, Vol. 114, No. 27, 2010
Wittich and Deiters
Figure 16. Dependence of the calculated Henry’s constant on the number of insertion cycles of carbon dioxide at 313 K in the (black line) orthorhombic case and the (red line) cubic simulation box. The relatively fast convergence against different limits is readily discerible.
TABLE 4: Resulting Energies of the Condensed Phase of [bmim][PF6] -1
T/K
∆VU/kJ mol
313.15 313.15 373.15 373.15
1822.8 ( 0.7 1822.5 ( 0.7 1816.6 ( 0.9 1817.5 ( 0.9
box geometry orthorh. cubic orthorh. cubic
box fluctuated around a cubic shape, but equilibration turned out to be much faster and the particle movement less hindered than with a cubic simulation box. Although properties like density and geometrical orientation can also be described by simulations with a cubic simulation box, sensitive energetic properties like Henry’s constant are described much better with a dynamic flexible box. This indicates that these properties are more influenced by finite-size effects and the boundary conditions. Furthermore, the diffusion of particles is hindered in the case of a cubic simulation box. This effect can be reproduced with different ensemble sizes. Of course, using a very large number of particles would decrease the influence of the simulation box, and at a certain system size it should become negligible, but this would use an excessive amount of computational time. Models that are considering Coulomb interactions explicitly yield good results for Henry’s constant when an orthorhombic simulation box is used. Although simple Lennard-Jones models are a good choice to describe many thermodynamic properties of the pure species, they cannot produce such good results for sensitive properties like Henry’s constant, because they do not consider the significant electrostatic contributions in ionic liquids. Whenever good gas models (with partial charges) were used, superior results were obtained. Nonpolar models probably do not give good estimations because of neglected polarizability effects. It is very probable that with increasing system size these effects decrease and eventually vanish for very large systems. From the simulations of Shah and Maginn1 it is obvious that the system size must be much larger than 250 cation-anion pairs. Because of our limited computing time, however, it was not possible to investigate the size effects in more detail. The data obtained with a smaller system of 64 cation-anion pairs were in excellent agreement with the results for the larger systems (cubic and noncubic).
It would be interesting to expand the method of a dynamic flexible to a Gibbs-ensemble to compare solubilities at different pressures with experimental data and data obtained in earlier simulations with a cubic box. With this new method it would also be possible to examine Henry’s constants of mixtures of different ionic liquids. References and Notes (1) Shah, J. K.; Maginn, E. J. Fluid Phase Equilib. 2004, 222, 195– 203. (2) Urukova, I.; Vorholz, J.; Maurer, G. J. Phys. Chem. B 2005, 109, 12154–12159. (3) Urukova, I.; Vorholz, J.; Maurer, G. J. Phys. Chem. B 2006, 110, 18072. (4) Wittich, B.; Deiters, U. K. Mol. Sim. 2010, 36, 364–372. (5) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021–12024. (6) Potoff, J. J.; Siepmann, J. I. AIChE J. 2001, 47, 1676–1682. (7) Cracknell, R. F. Phys. Chem. Chem. Phys. 2001, 3, 2091–2097. (8) Miyano, Y. Fluid Phase Equilib. 1999, 158-160, 29–35. (9) Bohn, M.; Lustig, R.; Fischer, J. Fluid Phase Equilib. 1986, 25, 251–262. (10) Straub, J. E.; Karplus, M. Chem. Phys. 1991, 158, 221–248. (11) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269–6271. (12) Vrabec, J.; Stoll, J.; Hasse, H. J. Phys. Chem. B 2001, 105, 12126– 12133. (13) Shing, K. S.; Gubbins, K. E.; Lucas, K. Mol. Phys. 1988, 65, 1235– 1252. (14) Widom, B. J. Stat. Phys. 1978, 19, 563–574. (15) Shing, K. S.; Chung, S. T. J. Phys. Chem. 1987, 91, 1674–1681. (16) Gu, Z.; Brennecke, J. F. J. Chem. Eng. Data 2002, 47, 339–345. (17) Tokuda, H.; Hayamizu, K.; Ishii, K.; Susan, M. A. B. H.; Watanabe, M. J. Phys. Chem. B 2004, 108, 16593–16600. ´ .; Tuma, D.; Xia, J.; Maurer, G. J. Chem. (18) Perez-Salado Kamps, A Eng. Data 2003, 48, 746–749. (19) Camper, D.; Becker, C.; Koval, C.; Noble, R. Ind. Eng. Chem. Res. 2005, 44, 1928–1933. ´ .; Tuma, D.; Maurer, G. Fluid (20) Kumełan, J.; Perez-Salado Kamps, A Phase Equilib. 2005, 228, 207–211. ´ .; Urukova, I.; Tuma, D.; (21) Kumełan, J.; Perez-Salado Kamps, A Maurer, G. J. Chem. Thermodyn. 2005, 37, 595–602. ´ .; Urukova, I.; Tuma, D.; (22) Kumełan, J.; Perez-Salado Kamps, A Maurer, G. J. Chem. Thermodyn. 2007, 39, 335. ´ .; Tuma, D.; Maurer, G. (23) Kumełan, J.; Perez-Salado Kamps, A J. Chem. Eng. Data 2006, 51, 11–14. (24) Anthony, J. L.; Maginn, E. J.; Brennecke, J. F. J. Phys. Chem. B 2002, 106, 7315–7320.
JP101676M
