Energetic Contributions from the Cation and Anion to the Stability of

Article pubs.acs.org/JPCB

Energetic Contributions from the Cation and Anion to the Stability of Carbon Dioxide Dissolved in Imidazolium-Based Ionic Liquids Ryosuke Ishizuka,†,‡ Nobuyuki Matubayasi,*,†,‡ Kai-Min Tu,§ and Yasuhiro Umebayashi∥ †

Division of Chemical Engineering, Graduate School of Engineering Science, Osaka University, Toyonaka, Osaka 560-8531, Japan Elements Strategy Initiative for Catalysts and Batteries, Kyoto University, Katsura, Kyoto 615-8520, Japan § Department of Chemistry, Graduate School of Science, Kyoto University, Kitashirakawa-Oiwakecho, Kyoto 606-8502, Japan ∥ Graduate School of Science and Technology, Niigata University, 8050, Ikarashi 2-no-cho, Nisi-ku, Niigata 950-2181, Japan ‡

ABSTRACT: Cation and anion effects were investigated in the energetics of CO2 solvation in room-temperature ionic liquids. The solvation free energy in 1-n-butyl-3methylimidazolium bis(trifluoromethylsulfonyl)imide ([C4mim][NTf2]) was calculated with three types of force fields by molecular dynamics simulations combined with the energy-representation (ER) method, and the interaction responsible for the CO2 stability was examined in terms of the average sum of the solute−solvent interaction energy and its electrostatic and van der Waals components. A key role of the van der Waals component was found for the cation contribution, and the electrostatic component was seen to be minor in the anion contribution. The solvation free energy of CO2 was also investigated in ionic liquids of the form 1-n-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([Cnmim][NTf2]) with varying alkyl-chain lengths of n = 4, 6, 8, and 12 by the ER method. The free energy was found to depend weakly on n, in agreement with experiments, and the cation and anion contributions and the electrostatic and van der Waals interaction components acted similarly at all values of n examined. The alkyl-chain length was found not to affect the local structure around CO2 strongly, and its effect on the interaction energy with CO2 appeared mainly through the bulk density.

■

INTRODUCTION Room-temperature ionic liquids (RTILs) have been paid much attention and regarded as novel solvents that can facilitate the separation of CO2 from gas mixtures.1−4 Many research groups have reported the Henry constants of CO2 in a variety of RTILs; this constant characterizes the CO2 solubility and is determined by the free energy of CO2 solvation in RTILs.5−20 A key to the rational design of an RTIL is the choice of cation and anion. Indeed, the dependence of the thermodynamics of CO2 solvation on the combination of cation and anion has been experimentally assessed, and a computational approach has also been employed to understand the intermolecular interactions of CO2 with RTILs.10,21−30 In general, RTILs are mixed solvents composed of large asymmetric organic cations and inorganic or organic anions. It is thus natural and often useful to separate the effects of the respective solvent species. The purpose of the present work is to examine the separate effects of the cation and anion on the energetics of CO2 solvation in RTILs. As noted in the preceding paragraph, the Henry constant is determined by the free energy of solvation, which reflects the local solvation structures in solution. Within the context of statistical thermodynamics, the ensemble of solution structures is integrated to the distribution of the solute−solvent interaction energies and further to the solvation free energy. The solute−solvent energetics is thus a useful measure in discussing solvation, particularly when the molecular structure is complex for the solute and/or solvent, © XXXX American Chemical Society

for example, for RTILs. In this work, we approach the energetics of CO2 solvation in RTILs using molecular dynamics (MD) simulations and statistical-mechanical method, by focusing on (1) 1-n-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([C4mim][NTf2]) with several force fields and (2) 1-n-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([Cnmim][NTf2]) with varying alkyl-chain lengths of n = 4, 6, 8, and 12. MD simulation is a powerful means for gaining molecularlevel information about intermolecular interactions between CO2 and RTILs. Analyses of thermodynamic quantities are possible with MD in terms of the separate contributions from cations and anions and further from such interaction components as electrostatic and van der Waals interactions. In the present work, these energetic contributions are discussed for CO2 solvation in RTILs in relation to the solvation free energy Δμtotal; Δμtotal is a thermodynamic observable governing the solubility. To determine Δμtotal, we employed the energy-representation (ER) theory, a density-functional method using distribution functions of solute−solvent interaction energy. In general, the separate contributions of cations and anions are not observable, and there is a need for a model of solvation to Received: October 9, 2014 Revised: December 29, 2014

A

DOI: 10.1021/jp5101957 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 1. Force Fields Employed Lennard-Jones parameter partial charge combining rule electrostatic 1−4 scaling factor van der Waals 1−4 scaling factor bonded interaction

model A

model B

model C

Köddermann et al.63 Lopes et al.64,65 geometric mean 1.0 0.5 Köddermann et al.63

Köddermann et al.63 Lopes et al.64,65 geometric mean 0.5 0.5 GAFF66,67

Lopes et al.64,65 Lopes et al.64,65 Lorentz−Berthelot 0.5 0.5 GAFF66,67

Ei = Eiel + Eivdw

conduct the separation. The ER method provides a straightforward way to decompose the solvation free energy Δμtotal into the separate contributions from cations and anions through31−35 Δμtotal = Δμcation + Δμanion

When the solute is at infinite dilution and the solute−solvent interaction is of site−site form, each component can be related to the local structure of the solvent around the solute by57 qαqj Eiel = 4πρ ∑ ∑ gαj(rαj)rαj 2 drαj r αj α ∈ solute j ∈ i (5)

(1)

∫

where Δμcation and Δμanion are the contributions from cations and anions, respectively. Widely used methods for computing the solvation free energy are the thermodynamic integration (TI) and free energy perturbation (FEP) methods. Although the TI and FEP methods are theoretically exact, numerical calculations are quite demanding to obtain the converged result. In contrast, the ER method reduces the computational cost significantly by formulating density-functional theory over the energy coordinate and employing an approximate functional for the solvation free energy. Among a variety of approximate methods for free energy,36−49 the ER method is unique in achieving a compromise among accuracy, efficiency, and range of applicability.32−34,50−56 However, the ER method has mainly been applied to aqueous systems. In the present work, we apply the ER method to CO2 in several imidazoliumbased ionic liquids and show the effectiveness of the use of the ER method for RTILs.

Eivdw = 4πρ

⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σαj σαj 4ϵαj⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥gαj(rαj)rαj 2 drαj r ⎝ rαj ⎠ ⎦ ⎣ ⎝ αj ⎠ (6)

where ρ is the number density of the solvent, rαj is the distance from the αth atom of the solute to the jth atom of the ith solvent species, gαj is the solute−solvent radial distribution function (RDF) between the αth and jth atoms, q is the atomic partial charge, ϵ and σ are Lennard-Jones parameters, and the sum is taken over all of the atoms within the solute and within the ith solvent species. Although E and its components are not observable, they are additive and thus straightforward for analyzing and interpreting the interaction components responsible for the stability of the solute in solution.34 Computational Procedures. We first focus on the energetics of CO2 solvation in [C4mim][NTf2] by examining three types of force fields. Because the force field for RTILs is still at the stage of active development,24,25,58−62 it is desirable to validate the MD results over different force fields. In this work, models A−C in Table 1 were employed for [C4mim][NTf2]. Model A is the force field developed by Köddermann et al.63 They refined the intramolecular and van der Waals parameters by Lopes et al.64,65 to reproduce experimental measurements of imidazolium-based ionic liquids.63 The atomic partial charges of [C4mim][NTf2] are the same as the original ones developed by Lopes et al.64,65 To investigate the influence of the bonded interaction parameters, we formulated model B by combining the nonbonded interaction parameters employed in model A with the bonded interaction parameters of the generalized Amber force field (GAFF) assigned by antechamber.66,67 Furthermore, to investigate the effects of the van der Waals parameters, model C employs the same parameters as in model B for the bonded interaction terms and the original parameters by Lopes et al.64,65 for the nonbonded terms. In a later part of the present article, the effect of the alkyl-chain length of the cation is studied, and model B was employed for [Cnmim][NTf2] with n = 6, 8, and 12. For all calculations, the force field of CO2 was taken from Harris and Yung.68 All of the simulations reported here were carried out at a temperature of 300 K and a pressure of 1 bar using the GROMACS 4.5.5 simulation program.69 Each simulation contained 216 ion pairs and a single CO2 molecule in a cubic unit cell with periodic boundary conditions and was conducted over 100 ns with the leapfrog stochastic dynamics integrator at

METHODS Theoretical Background. In the energy-representation (ER) method, the contribution Δμi (where i denotes either cation or anion) from each solvent species is written for the pairwise-additive solute−solvent interaction as

∫ dϵi f (ϵi) = ∫ dϵi ϵiρi (ϵi) + ∫ dϵi f (ϵi) (2)

where ϵi is the pair interaction energy between the solute and the ith solvent species, Ei is the average sum of the interaction energy with the ith solvent species in the solution system of interest, ρi(ϵi) is the average distribution (histogram) of the pair interaction energy ϵi in the solution, and f(ϵi) is a functional taking into account the effect of solvent reorganization including the excluded-volume effect.34,35 In the present work, the approximate functional described in previous works32,35,52 was employed; actually, only the second term of eq 2 is treated approximately, and the first term is evaluated exactly under a given set of force fields. The average interaction energy Etotal is the sum of the cation and anion contributions and is expressed as Etotal = Ecation + Eanion

∑ ∑∫ α ∈ solute j ∈ i

■

Δμi = Ei +

(4)

(3)

When the solute−solvent interaction consists of electrostatic and van der Waals potentials, Ei is further decomposed into the electrostatic component Eeli and the van der Waals component Evdw through i B


Article

The Journal of Physical Chemistry B

According to Figure 1, the solvation free energy Δμtotal and the contributions from the cation and anion are quantitatively similar between models A and B, showing an insensitivity to the examined choice of bonded interactions. The comparison between models B and C indicates that the van der Waals (Lennard-Jones) interactions between CO2 and RTILs affect the CO2 stabilization more strongly. Actually, we will see later that the van der Waals interactions play a key role in the energetics of CO2 solvation. Figure 1 further shows a plot of the average sum of the solute−solvent interaction energy Etotal and its cation and anion contributions. It is evident that Δμtotal and Etotal are more negative for model C than for models A and B. The cation contributions are seen to be similar among all the models, and the anion contributions are responsible for the model dependence. On the basis of the parallelism of the model dependence of Δμ and E, we focus on E and its components. As noted with respect to eqs 2−4, E is additive and is more straightforward than Δμ for discussing the interaction component responsible for the stability of CO2 in RTILs. In addition, the cation has a positively charged imidazolium ring and a nearly neutral alkyl chain. It will then be insightful to decompose Ecation further into the contributions from the ring Eimid and from the alkyl chain Ealkyl and to examine the separate effects of the electrostatic and van der Waals interaction components through eq 4. The definitions of the imidazolium-ring and alkyl-chain groups are drawn in Figure 2, and the solute−solvent energy E

a time step of 2 fs and an inverse friction constant of 0.2 ps.70 The pressure was controlled using the Parrinello−Rahman barostat at a coupling time of 1 ps. The electrostatic interactions were computed using the smooth particle mesh Ewald (PME) method71 with a real-space cutoff of 13.5 Å; a spline order of 6; a relative tolerance of 10−5 (inverse decay length of 0.23 Å−1); and a reciprocal-space mesh size of 32 for each of the x, y, and z directions. The van der Waals interactions were cut off beyond 12.0 Å, with switching range of 10.0−12.0 Å. The long-range correction of the LJ interaction was applied after the MD simulation to the solvation free energy and the average sum of the solute−solvent interaction energy (solute−solvent energy) using the standard scheme72,73 by supposing that the radial distribution functions are unity beyond the switching region. All bond lengths with hydrogen atoms were held fixed using the LINCS algorithm.74 The solvation free energy was calculated by the energy-representation (ER) method, and the method also requires the MD trajectories of the pure solvent and the isolated solute. The simulation conditions for pure RTIL were the same as those for the solution MD, except that the MD length was 10 ns, and the single-molecule simulation for isolated CO2 was carried out in a vacuum for 100 ns. The configuration was sampled every 1 ps for pure RTIL and every 100 fs for isolated CO2. The testparticle insertion of CO2 was carried out 100 times per puresolvent configuration sampled, leading to the generation of 106 solute−solvent configurations in total. See previous articles for methodological details and the explicit form of the functional for the solvation free energy.32,35,52 The first term of eq 2 is the solute−solvent energy, and its components are described by eqs 3−6 in the Theoretical Background subsection. They are computed only from the MD simulations of the solution system of interest, without any reference to the test-particle insertion, which employs the trajectories of the pure solvent and the isolated solute.57

■

Figure 2. Ecation decomposed into contributions from (a) the imidazolium ring and (b) the alkyl chain. The methyl and methylene groups bonded to the ring are included in the imidazolium-ring group. The other hydrocarbon parts are counted as the alkyl-chain group.

RESULTS AND DISCUSSION Force-Field Dependence. In Figure 1, Δμtotal and the contributions from the solvent species in [C4mim][NTf2] are

and its electrostatic and van der Waals components are shown in Figure 3. For all of the models, Eelanion is relatively small in magnitude, and Evdw anion is the major component of Eanion. To connect the energetics of CO2−anion interactions with the local structure of solvation, in Figure 4 we show the radial distribution functions (RDFs) averaged over the heavy (nonhydrogen) atoms in the imidazolium ring, alkyl chain, and anion with respect to the carbon and oxygen atoms of CO2.

Figure 1. Δμtotal and Etotal calculated for models A−C with the cation and anion contributions. The error bars are shown at the 95% confidence level.70

plotted against the employed models. The Δμtotal value from ER (−1.5 kcal/mol) with model A is in agreement with the experimental (−0.6 kcal/mol) and computational (−0.7 kcal/ mol) results;10,27 model A was built to fit experimental quantities. The approximate functional employed in the ER formalism can thus be useful for the study of CO2 solvation in RTILs.70

Figure 3. Electrostatic and van der Waals components of Ei (i = imidazolium ring, alkyl chain, and anion) for models A−C. The error bars are shown at the 95% confidence level.70 C


Article


Figure 4. Radial distribution functions (RDFs) obtained for models A (blue), B (red), and C (green). The top, center, and bottom plots refer to the imidazolium ring of [C4mim], the alkyl chain of [C4mim], and [NTf2], respectively, and are the averages of the RDFs between the heavy (nonhydrogen) atoms in the respective groups and the carbon (left column) and oxygen (right column) atoms of CO2.

leading to a favorable (negative) Eelimid value. As for the alkyl chain, the electrostatic component Eelalkyl in Figure 3 is close to zero, as expected from the small charges assigned to the hydrocarbon chain. The RDFs in Figure 4c,d have distinct peaks, and the van der Waals interaction with the alkyl chain contributes to the stabilization of CO2. In summary, the electrostatic component is effective only for the interaction of the solute CO2 with the imidazolium ring, and the van der Waals component Evdw total in the total solute−solvent energy is a few times larger than the electrostatic component Eeltotal. This shows the key role of van der Waals interactions for CO2 stabilization in imidazolium-based ionic liquids, in agreement with refs 26 and 27. Effect of Alkyl-Chain Length. We now turn to the effect of the alkyl-chain length in the imidazolium-based cation. Using model B in Table 1, we examine the energetics of CO2 solvation in [Cnmim][NTf2] with n = 4, 6, 8, and 12. The densities of [Cnmim][NTf2] with model B are listed in Table 2, in good agreement with the experimental values.75 In Figure 5, the solvation free energy Δμtotal and its cation and anion contributions are shown with experimental data76 of Δμtotal. For the alkyl-chain lengths of n = 4, 6, and 8, the solvation free energies from the ER method are in agreement with the experimental values. It is further seen in both the theoretical and experimental results that the effect of the alkyl-chain length on the CO2 Δμ in imidazolium-based RTILs is weak.27 According to the decomposition into cation and anion contributions, the cation is responsible for the favorable

Table 2. Mass Densities of [Cnmim][NTf2] and Molar Concentrations of the Chemical Groups with n = 4, 6, 8, and 12 and the Corresponding Experimental Values75 mass density (g/cm3)

[C4mim] [NTf2] [C6mim] [NTf2] [C8mim] [NTf2] [C12mim] [NTf2]

molar concentration (mol/L)b

MD simulationa

experiment

imidazolium ring

alkyl chain

anion

1.435

1.436

23.97

10.27

51.37

1.368

1.371

21.42

15.30

45.91

1.317

1.319

19.41

19.41

41.59

1.241

1.245

16.36

25.71

35.06

Error at the 95% confidence level is less than ±0.001 g/cm3. bSum of the average number densities of the heavy (non-hydrogen) atoms involved in the respective group.

a

The RDFs of the anion are almost flat around both carbon and oxygen. Equation 5 then leads to the small Eelanion value in Figure 3 due to the charge neutrality of CO2; the repulsive electrostatic interaction between the oxygen atoms of CO2 and the anion reduces Eelanion (in magnitude). In contrast to Eanion, the electrostatic component Eelimid and the van der Waals component Evdw imid for the imidazolium ring contribute equally to Eimid. Although the RDF of the imidazolium ring around the carbon atom of CO2 in Figure 4a has a distinct peak, Figure 4b shows that the negatively charged oxygen atom comes closer to the imidazolium ring, D


Article


Figure 5. Δμtotal and Etotal for [Cnmim][NTf2] with n = 4, 6, 8, and 12 with the cation and anion contributions. The error bars are shown at the 95% confidence level.70 The data for n = 4 are the same as those in Figure 1 for model B. The experimental data for the solvation free energy Δμexp were taken from ref 76 and were absent for n = 12 to the best of our knowledge.

(negative) Δμtotal value, and the anion contribution is marginally unfavorable (positive). The dependence on the alkyl-chain length is also weak for the cation and anion contributions. Figure 5 also shows the average sum of the solute−solvent interaction energy Etotal and its cation and anion contributions. The parallelism is then observed between Δμ and E over the range of n values examined, and we focus on E and its components in the following discussion. Figure 6 provides the

(relatively) strong variation with n corresponds to the change in the molar concentration in Table 2. According to Figure 6, the decrease of Evdw alkyl counteracts the el increases of Evdw imid and Eimid. This compensation is partial, and Ecation of Figure 5 varies with the alkyl-chain length n by several tenths of a kilocalorie per mole. The n dependence is then 27 canceled between Ecation and Evdw anion, leading to the insensitivity of Etotal to n in Figure 5, which corresponds, in turn, to the invariable Δμtotal. Our final remark is about the effect of excluded volume. This effect is the major part of the repulsive interaction between solute and solvent and is the cause of the free-energy penalty for removing the solvent from the region into which the solute is to be inserted. In the ER method, the excluded-volume component Δμexcl can be obtained by restricting the integral over ϵi in eq 2 to a high-energy domain ranging from ϵci to infinity, where ϵci is a threshold value for defining the excludedvolume domain. The threshold value is taken to satisfy a requirement that the domain of ϵi > ϵci corresponds to the solute−solvent overlap and is essentially inaccessible in the solution system of interest. We set ϵci = 10 kcal/mol, although the following discussion is valid irrespective of the (reasonable) setting of the ϵci value. It should be noted that Δμexcl is a useful but nonobservable quantity. A model of solvation needs to be employed to introduce Δμexcl. In the present work, Δμexcl is discussed within the ER formalism, and its cation and anion contributions can be examined by virtue of eqs 1 and 2. The excluded-volume component Δμexcl is presented in Figure 8 with its cation and anion contributions. It is observed that Δμexcl total stays constant with the variation of the alkyl-chain length n. Etotal given in Figure 5 captures the effect of attractive excl interactions between the solute and solvent, and Δμtotal represents the repulsive-interaction effect. Figures 5 and 8 show that both depend weakly on n. The weak dependence of Δμexcl total is due to the compensating dependence on n of the excl cation contribution Δμexcl cation and the anion contribution Δμanion, as was also the case for the solute−solvent energy E discussed excl above with respect to Figures 5 and 6. Δμexcl cation and Δμanion vary by a few tenths of a kilocalorie per mole in opposite directions to Ecation and Eanion in Figure 5, furthermore leading to the cancellation of the n dependence between Ei and Δμexcl (i = i cation or anion) and to the weak effect of n on Δμi.

Figure 6. Electrostatic and van der Waals components of E for the imidazolium ring, alkyl chain, and anion in [Cnmim][NTf2] with n = 4, 6, 8, and 12. The error bars are shown at the 95% confidence level.70 The data for n = 4 are the same as those in Figure 3 for model B.

electrostatic and van der Waals components of Eimid, Ealkyl, and Eanion at n = 4, 6, 8, and 12, where the imidazolium-ring and alkyl-chain groups are defined in Figure 2. For the following discussion of the local structures of the cation and anion around CO2, Figure 7 also shows plots of the RDFs of the imidazolium ring, alkyl chain, and anion. The energetics presented in Figure 6 with varying length of alkyl chain behaves similarly to that shown in Figure 3 for n = 4. Eelanion is close to zero, and Evdw anion is the dominant component in el Eanion. For the cation, Evdw imid and Eimid contribute equally to Eimid, el and Ealkyl is small, as expected from the charges on the el vdw hydrocarbon. The n dependence is evident for Evdw imid, Eimid, Ealkyl, vdw and Eanion. This is not due to the response of the local structure to the chain length n, however. The RDFs for the imidazolium ring and anion do not change with n, and eqs 5 and 6 then el vdw show that Evdw imid, Eimid, and Eanion vary with n in correspondence to the molar concentration for the (heavy) atoms in the respective groups listed in Table 2. The n dependence is observed in Figure 7 for the RDFs of the alkyl-chain group; the first peak is less distinct when the alkyl chain is longer. Still, Evdw alkyl is affected more strongly by the bulk density of the methylene and methyl groups forming the alkyl chain, and its

■

CONCLUSIONS The energetics of CO2 solvation in imidazolium-based ionic liquids [Cnmim][NTf2] were investigated over a wide range of alkyl-chain lengths n = 4, 6, 8, and 12. The solvation free energy was computed by the energy-representation method, in E


Article


Figure 7. Radial distribution functions (RDFs) for n = 4 (red), 6 (blue), 8 (yellow), and 12 (green). The top, center, and bottom plots refer to the imidazolium ring of [Cnmim], the alkyl chain of [Cnmim], and [NTf2], respectively, and are the averages of the RDFs between the heavy (nonhydrogen) atoms in the respective groups and the carbon (left column) and oxygen (right column) atoms of CO2. The data for n = 4 are the same as those in Figure 4 for model B.

of CO2 mainly through the bulk density. In our analysis, the separate effects of cation and anion were treated and the role of van der Waals interactions was highlighted. The analysis scheme is not restricted to CO2 solvation in [Cnmim][NTf2] and will be employed in subsequent work for a general class of gas solvation in RTILs.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

Figure 8. Excluded-volume component of the solvation free energy with the cation and anion components.

Notes

The authors declare no competing financial interest.

■

agreement with experiments over the n values examined. The decomposition into the cation and anion contributions was then conducted for the solvation free energy and the average sum of the solute−solvent interaction energy. It was found that the anion contribution governs the dependence on the force field and that the solvation free energy and solute−solvent energy behave in parallel for the dependence on the alkyl-chain length. Within the solute−solvent energy, the van der Waals component plays the key role. The electrostatic component is minor for the anion. The length of the alkyl chain attached to the imidazolium ring does not have a strong effect on the structural features observed for RDFs and affects the interaction

ACKNOWLEDGMENTS This work was supported by Grants-in-Aid for Scientific Research (Nos. 23651202 and 26240045) from the Japan Society for the Promotion of Science; by the Elements Strategy Initiative for Catalysts and Batteries from the Ministry of Education, Culture, Sports, Science, and Technology; and by the Computational Materials Science Initiative, Theoretical and Computational Chemistry Initiative, and Strategic Programs for Innovative Research of the Next-Generation Supercomputing Project. The simulations were conducted partly using computational resources of the HPCI systems provided by COMA at F


Article


(20) Gonzalez-Miquel, M.; Bedia, J.; Abrusci, C.; Palomar, J.; Rodriguez, F. Anion Effects on Kinetics and Thermodynamics of CO2 Absorption in Ionic Liquids. J. Phys. Chem. B 2013, 117, 3398−3406. (21) Kanakubo, M.; Umecky, T.; Hiejima, Y.; Aizawa, T.; Nanjo, H.; Kameda, Y. Solution Structures of 1-Butyl-3-methylimidazolium Hexafluorophosphate Ionic Liquid Saturated with CO2: Experimental Evidence of Specific Anion−CO2 Interaction. J. Phys. Chem. B 2005, 109, 13847−13850. (22) Shah, J. K.; Maginn, E. J. Monte Carlo Simulations of Gas Solubility in the Ionic Liquid 1-n-Butyl-3-methylimidazolium Hexafluorophosphate. J. Phys. Chem. B 2005, 109, 10395−10405. (23) Huang, X.; Margulis, C. J.; Li, Y.; Berne, B. J. Why Is the Partial Molar Volume of CO2 So Small When Dissolved in a Room Temperature Ionic Liquid? Structure and Dynamics of CO2 Dissolved in [Bmim+] [PF6−]. J. Am. Chem. Soc. 2005, 127, 17842−17851. (24) Bhargava, B. L.; Balasubramanian, S. Probing Anion−Carbon Dioxide Interactions in Room Temperature Ionic Liquids: Gas Phase Cluster Calculations. Chem. Phys. Lett. 2007, 444, 242−246. (25) Bhargava, B. L.; Balasubramanian, S. Insights into the Structure and Dynamics of a Room-Temperature Ionic Liquid: Ab Initio Molecular Dynamics Simulation Studies of 1-n-Butyl-3-methylimidazolium Hexafluorophosphate ([bmim][PF6]) and the [bmim][PF6]− CO2 Mixture. J. Phys. Chem. B 2007, 111, 4477−4487. (26) Shi, W.; Maginn, E. J. Molecular Simulation and Regular Solution Theory Modeling of Pure and Mixed Gas Absorption in the Ionic Liquid 1-n-Hexyl-3-methylimidazolium Bis(trifluoromethylsulfonyl)amide ([hmim][Tf2N]). J. Phys. Chem. B 2008, 112, 16710−16720. (27) Kerlé, D.; Ludwig, R.; Geiger, A.; Paschek, D. Temperature Dependence of the Solubility of Carbon Dioxide in Imidazolium-Based Ionic Liquids. J. Phys. Chem. B 2009, 113, 12727−12735. (28) Sambasivarao, S. V.; Acevedo, O. Development of OPLS-AA Force Field Parameters for 68 Unique Ionic Liquids. J. Chem. Theory Comput. 2009, 5, 1038−1050. (29) Kanakubo, M.; Aizawa, T.; Nanjo, H.; Kameda, Y.; Amo, Y.; Usuki, T. Liquid Structures of 1-Butyl-3-methylimidazolium Tetrafluoroborate and Carbon Dioxide Mixtures by X-ray Diffraction Measurements. Fluid Phase Equilib. 2010, 297, 183−186. (30) Lourenco, T. C.; Coelho, M. F. C.; Ramalho, T. C.; van der Spoel, D.; Costa, L. T. Insights on the Solubility of CO2 in 1-Ethyl-3methylimidazolium Bis(trifluoromethylsulfonyl)imide from the Microscopic Point of View. Environ. Sci. Technol. 2013, 47, 7421−7429. (31) Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energetic Representation. I. Formulation. J. Chem. Phys. 2000, 113, 6070−6081. (32) Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energy Representation. II. Functional for the Chemical Potential. J. Chem. Phys. 2002, 117, 3605−3616; J. Chem. Phys. 2003, 118, 2446 (erratum).. (33) Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energy Representation. III. Treatment of the Molecular Flexibility. J. Chem. Phys. 2003, 119, 9686−9702. (34) Karino, Y.; Matubayasi, N. Interaction-Component Analysis of the Urea Effect on Amino Acid Analogs. Phys. Chem. Chem. Phys. 2013, 15, 4377−4391. (35) Sakuraba, S.; Matubayasi, N. ERmod: Fast and Versatile Computation Software for Solvation Free Energy with Approximate Theory of Solutions. J. Comput. Chem. 2014, 35, 1592−1608. (36) Levy, R. M.; Belhadj, M.; Kitchen, D. B. Gaussian Fluctuation Formula for Electrostatic Free-Energy Changes in Solution. J. Chem. Phys. 1991, 95, 3627−3633. (37) Luzhkov, V.; Warshel, A. Microscopic Models for Quantum Mechanical Calculations of Chemical Processes in Solutions: LD/ AMPAC and SCAAS/AMPAC Calculations of Solvation Energies. J. Comput. Chem. 1992, 13, 199−213. (38) Åqvist, J.; Medina, C.; Samuelsson, J.-E. A New Method for Predicting Binding Affinity in Computer-Aided Drug Design. Protein Eng. 1994, 7, 385−391.

University of Tsukuba, TSUBAME2.5 at Tokyo Institute of Technology, and Cray XC30 at Kyoto University through the HPCI System Research Project (Projects hp120093, hp130022, hp140156, and hp140214). K.-M.T. is also grateful for a fellowship from the Interchange Association, Tokyo, Japan.

■

REFERENCES

(1) Anthony, J. L.; Maginn, E. J.; Brennecke, J. F. Solubilities and Thermodynamic Properties of Gases in the Ionic Liquid 1-n-Butyl-3methylimidazolium Hexafluorophosphate. J. Phys. Chem. B 2002, 106, 7315−7320. (2) Rogers, R. D.; Seddon, K. R. Ionic LiquidsSolvents of the Future? Science 2003, 302, 792−793. (3) Endres, F.; Zein El Abedin, S. Air and Water Stable Ionic Liquids in Physical Chemistry. Phys. Chem. Chem. Phys. 2006, 8, 2101−2116. (4) Zhang, X.; Zhang, X.; Dong, H.; Zhao, Z.; Zhang, S.; Huang, Y. Carbon Capture with Ionic Liquids: Overview and Progress. Energy Environ. Sci. 2012, 5, 6668−6681. (5) Blanchard, L. A.; Gu, Z.; Brennecke, J. F. High-Pressure Phase Behavior of Ionic Liquid/CO2 Systems. J. Phys. Chem. B 2001, 105, 2437−2444. (6) Cammarata, L.; Kazarian, S. G.; Salter, P. A.; Welton, T. Molecular States of Water in Room Temperature Ionic Liquids. Phys. Chem. Chem. Phys. 2001, 3, 5192−5200. (7) Cadena, C.; Anthony, J. L.; Shah, J. K.; Morrow, T. I.; Brennecke, J. F.; Maginn, E. J. Why Is CO2 So Soluble in Imidazolium-Based Ionic Liquids? J. Am. Chem. Soc. 2004, 126, 5300−5308. (8) Aki, S. N. V. K.; Mellein, B. R.; Saurer, E. M.; Brennecke, J. F. High-Pressure Phase Behavior of Carbon Dioxide with ImidazoliumBased Ionic Liquids. J. Phys. Chem. B 2004, 108, 20355−20365. (9) Shariati, A.; Peters, C. J. High-Pressure Phase Behavior of Systems with Ionic Liquids: II. The Binary System Carbon Dioxide + 1-Ethyl-3-methylimidazolium Hexafluorophosphate. J. Supercrit. Fluids 2004, 29, 43−48. (10) Anthony, J. L.; Anderson, J. L.; Maginn, E. J.; Brennecke, J. F. Anion Effects on Gas Solubility in Ionic Liquids. J. Phys. Chem. B 2005, 109, 6366−6374. (11) Kim, Y. S.; Choi, W. Y.; Jang, J. H.; Yoo, K. P.; Lee, C. S. Solubility Measurement and Prediction of Carbon Dioxide in Ionic Liquids. Fluid Phase Equilib. 2005, 228, 439−445. (12) Shiflett, M. B.; Yokozeki, A. Solubilities and Diffusivities of Carbon Dioxide in Ionic Liquids: [bmim][PF6] and [bmim][BF4]. Ind. Eng. Chem. Res. 2005, 44, 4453−4464. (13) Jacquemin, J.; Husson, P.; Majer, V.; Gomes, M. F. C. Influence of the Cation on the Solubility of CO2 and H2 in Ionic Liquids Based on the Bis(trifluoromethylsulfonyl)imide Anion. J. Solution Chem. 2007, 36, 967−979. (14) Kim, Y. S.; Jang, J. H.; Lim, B. D.; Kang, J. W.; Lee, C. S. Solubility of Mixed Gases Containing Carbon Dioxide in Ionic Liquids: Measurements and Predictions. Fluid Phase Equilib. 2007, 256, 70−74. (15) Muldoon, M. J.; Aki, S. N. V. K.; Anderson, J. L.; Dixon, J. K.; Brennecke, J. F. Improving Carbon Dioxide Solubility in Ionic Liquids. J. Phys. Chem. B 2007, 111, 9001−9009. (16) Yuan, X.; Zhang, S.; Lid, J.; Lu, X. Solubilities of CO2 in Hydroxyl Ammonium Ionic Liquids at Elevated Pressures. Fluid Phase Equilib. 2007, 257, 195−200. (17) Kurnia, K.; Harris, F.; Wilfred, C.; Mutalib, M. A.; Murugesan, T. Thermodynamic Properties of CO2 Absorption in Hydroxyl Ammonium Ionic Liquids at Pressures of (100−1600) kPa. J. Chem. Thermodyn. 2009, 41, 1069−1073. (18) Bara, J. E.; Carlisle, T. K.; Gabriel, C. J.; Camper, D.; Finotello, A.; Gin, D. L.; Noble, R. D. Guide to CO2 Separations in ImidazoliumBased Room-Temperature Ionic Liquids. Ind. Eng. Chem. Res. 2009, 48, 2739−2751. (19) Supasitmongkol, S.; Styring, P. High CO2 Solubility in Ionic Liquids and a Tetraalkylammonium-based Poly(ionic liquid). Energy Environ. Sci. 2010, 3, 1961−1972. G


Article


preferable to determine Eeli by integrating the RDF multiplied by the charge−charge interaction over the space. In practice, Eiel was calculated by directly computing the electrostatic part of the solute− solvent interaction at each snapshot and averaging over the trajectory. (58) Liu, Z.; Huang, S.; Wang, W. A Refined Force Field for Molecular Simulation of Imidazolium-Based Ionic Liquids. J. Phys. Chem. B 2004, 108, 12978−12989. (59) Bhargava, B. L.; Balasubramanian, S.; Klein, M. L. Modelling Room Temperature Ionic Liquids. Chem. Commun. 2008, 3339−3351. (60) Bhargava, B. L.; Saharay, M.; Balasubramanian, S. Ab Initio Studies on [bmim][PF6]−CO2 Mixture and CO2 Clusters. Bull. Mater. Sci. 2008, 31, 327−334. (61) Schmidt, J.; Krekeler, C.; Dommert, F.; Zhao, Y.; Berger, R.; Site, L. D.; Holm, C. Ionic Charge Reduction and Atomic Partial Charges from First-Principles Calculations of 1,3-Dimethylimidazolium Chloride. J. Phys. Chem. B 2010, 114, 6150−6155. (62) Dommert, F.; Wendler, K.; Berger, R.; Delle Site, L.; Holm, C. Force Fields for Studying the Structure and Dynamics of Ionic Liquids: A Critical Review of Recent Developments. ChemPhysChem 2012, 13, 1625−1637. (63) Köddermann, T.; Paschek, D.; Ludwig, R. Molecular Dynamic Simulations of Ionic Liquids: A Reliable Description of Structure, Thermodynamics and Dynamics. ChemPhysChem 2007, 8, 2464− 2470. (64) Lopes, J. N. C.; Deschamps, J.; Pádua, A. A. H. Modeling Ionic Liquids Using a Systematic All-Atom Force Field. J. Phys. Chem. B 2004, 108, 2038−2047. (65) Lopes, J. N. C.; Pádua, A. A. H. Molecular Force Field for Ionic Liquids Composed of Triflate or Bistriflylimide Anions. J. Phys. Chem. B 2004, 108, 16893−16898. (66) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (67) Sousa da Silva, A.; Vranken, W. ACPYPE - AnteChamber PYthon Parser interfacE. BMC Res. Notes 2012, 5, 367. (68) Harris, J. G.; Yung, K. H. Carbon Dioxide’s Liquid-Vapor Coexistence Curve And Critical Properties as Predicted by a Simple Molecular Model. J. Phys. Chem. 1995, 99, 12021−12024. (69) 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, 435−447. (70) The long MD simulation is needed to obtain converged radial distribution functions. The solvation free energy Δμtotal and its cation and anion components converge to errors at the 95% confidence level of 0.1 kcal/mol within 10 ns, except for the anion component Δμanion in Figure 1 for model C, which has an error of 0.2 kcal/mol. When the solution MD simulation is extended to 100 ns, the errors at the 95% confidence level of Δμtotal, Δμcation, and Δμanion are smaller than 0.1 kcal/mol. The average sum of the solute−solvent interaction energy (solute−solvent energy) E after a 100-ns MD simulation of the solution system is mostly within an error at the 95% confidence level of 0.1 kcal/mol, with 0.2−0.3 kcal/mol for some cases. In the method of energy representation, the solvation free energy converges faster than the average sum of the solute−solvent interaction energy, as noted in refs 35, 52, and 54. This is due to the variational principle, and the convergence of Δμtotal is faster than those of Δμcation, Δμanion, Etotal, Ecation, Eanion, and their electrostatic and van der Waals components. (71) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (72) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (73) Shirts, M. R.; Pitera, J. W.; Swope, W. C.; Pande, V. S. Extremely Precise Free Energy Calculations of Amino Acid Side Chain Analogs: Comparison of Common Molecular Mechanics Force Fields For Proteins. J. Chem. Phys. 2003, 119, 5740−5761.

(39) Carlson, H. A.; Jorgensen, W. L. An Extended Linear Response Method for Determining Free Energies of Hydration. J. Phys. Chem. 1995, 99, 10667−10673. (40) Kast, S. M. Combinations of Simulation and Integral Equation Theory. Phys. Chem. Chem. Phys. 2001, 3, 5087−5092. (41) Vener, M. V.; Leontyev, I. V.; Dyakov, Y. A.; Basilevsky, M. V.; Newton, M. D. Application of the Linearized MD Approach for Computing Equilibrium Solvation Free Energies of Charged and Dipolar Solutes in Polar Solvents. J. Phys. Chem. B 2002, 106, 13078− 13088. (42) Hirata, F., Ed. Molecular Theory of Solvation; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2003. (43) Galván, I. F.; Sánchez, M. L.; Martín, M. E.; Olivares del Valle, F. J.; Aguilar, M. A. Geometry Optimization of Molecules in Solution: Joint Use of the Mean Field Approximation and the Free-Energy Gradient Method. J. Chem. Phys. 2003, 118, 255−263. (44) Freedman, H.; Truong, T. N. Coupled Reference Interaction Site Model/Simulation Approach for Thermochemistry of Solvation: Theory and Prospects. J. Chem. Phys. 2004, 121, 2187−2198. (45) Higashi, M.; Hayashi, S.; Kato, S. Geometry Optimization Based on Linear Response Free Energy with Quantum Mechanical/ Molecular Mechanical Method: Applications to Menshutkin-type and Claisen Rearrangement Reactions in Aqueous Solution. J. Chem. Phys. 2007, 126, 144503. (46) Chuev, G. N.; Fedorov, M. V.; Crain, J. Improved Estimates for Hydration Free Energy Obtained by the Reference Interaction Site Model. Chem. Phys. Lett. 2007, 448, 198−202. (47) Yamamoto, T. Variational and Perturbative Formulations of Quantum Mechanical/Molecular Mechanical Free Energy with MeanField Embedding and Its Analytical Gradients. J. Chem. Phys. 2008, 129, 244104. (48) Frolov, A. I.; Ratkova, E. L.; Palmer, D. S.; Fedorov, M. V. Hydration Thermodynamics Using the Reference Interaction Site Model: Speed or Accuracy? J. Phys. Chem. B 2011, 115, 6011−6022. (49) Kokubo, H.; Hu, C. Y.; Pettitt, B. M. Peptide Conformational Preferences in Osmolyte Solutions: Transfer Free Energies of Decaalanine. J. Am. Chem. Soc. 2011, 133, 1849−1858. (50) Takahashi, H.; Matubayasi, N.; Nakahara, M.; Nitta, T. A Quantum Chemical Approach to the Free Energy Calculations in Condensed Systems: The QM/MM Method Combined with the Theory of Energy Representation. J. Chem. Phys. 2004, 121, 3989− 3999. (51) Matubayasi, N.; Liang, K. K.; Nakahara, M. Free-Energy Analysis of Solubilization in Micelle. J. Chem. Phys. 2006, 124, 154908. (52) Matubayasi, N.; Shinoda, W.; Nakahara, M. Free-Energy Analysis of the Molecular Binding into Lipid Membrane with the Method of Energy Representation. J. Chem. Phys. 2008, 128, 195107. (53) Takahashi, H.; Ohno, H.; Kishi, R.; Nakano, M.; Matubayasi, N. Computation of the Free Energy Change Associated with OneElectron Reduction of Coenzyme Immersed in Water: A Novel Approach within the Framework of the Quantum Mechanical/ Molecular Mechanical Method Combined with the Theory of Energy Representation. J. Chem. Phys. 2008, 129, 205103. (54) Karino, Y.; Fedorov, M. V.; Matubayasi, N. End-Point Calculation of Solvation Free Energy of Amino-Acid Analogs by Molecular Theories of Solution. Chem. Phys. Lett. 2010, 496, 351−355. (55) Kawakami, T.; Shigemoto, I.; Matubayasi, N. Free-Energy Analysis of Water Affinity in Polymer Studied by Atomistic Molecular Simulation Combined with the Theory of Solutions in the Energy Representation. J. Chem. Phys. 2012, 137, 234903. (56) Takemura, K.; Guo, H.; Sakuraba, S.; Matubayasi, N.; Kitao, A. Evaluation of Protein−Protein Docking Model Structures Using AllAtom Molecular Dynamics Simulations Combined with the Solution Theory in the Energy Representation. J. Chem. Phys. 2012, 137, 215105. (57) The computation of Eeli is not done with the Coulomb interaction in its bare form of 1/r, but with the smooth PME form described in ref 35, because our MD approach employs the smooth PME to handle the electrostatic interaction. Numerically, it is not H


Article

The Journal of Physical Chemistry B (74) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (75) Tariq, M.; Forte, P. A. S.; Gomes, M. F. C.; Lopes, J. N. C.; Rebelo, L. P. N. Densities and Refractive Indices of Imidazolium- and Phosphonium-based Ionic Liquids: Effect of Temperature, Alkyl Chain Length, and Anion. J. Chem. Thermodyn. 2009, 41, 790−798. (76) Baltus, R. E.; Culbertson, B. H.; Dai, S.; Luo, H.; DePaoli, D. W. Low-Pressure Solubility of Carbon Dioxide in Room-Temperature Ionic Liquids Measured with a Quartz Crystal Microbalance. J. Phys. Chem. B 2004, 108, 721−727.

I


Energetic Contributions from the Cation and Anion to the Stability of

Recommend Documents