Article pubs.acs.org/JPCB
Contribution of the Acetate Anion to CO2 Solubility in Ionic Liquids: Theoretical Method Development and Experimental Study Wei Shi,*,†,‡,§ Robert L. Thompson,†,‡ Erik Albenze,†,‡ Janice A. Steckel,† Hunaid B. Nulwala,∥,† and David R. Luebke† †
U.S. Department of Energy, National Energy Technology Laboratory, 626 Cochrans Mill Road, Pittsburgh, Pennsylvania 15236, United States ‡ URS Corporation, South Park, Pennsylvania 15129, United States § Department of Chemical and Petroleum Engineering, University of Pittsburgh, 3700 O’Hara Street, Pittsburgh, Pennsylvania 15261, United States ∥ Department of Chemistry, Carnegie Mellon University, 4400 Fifth Avenue, Pittsburgh, Pennsylvania 15213, United States ABSTRACT: A new theoretical method was developed to compute the Henry’s law constant for gas absorption in a solvent through strong nonphysical interactions. The new method was created by expanding the test particle insertion method typically applied to physisorbing systems to account for the strong intermolecular interactions present in chemisorbing systems. By using an ab initio (AI)-based Boltzmann-averaged potential to model the interaction between CO2 and the tetra-n-butylphosphonium acetate ([P4444][CH3COO]) ionic liquid, the total Henrys’s law constant at 298 K was computed to be 0.011 to 0.039 bar, reasonably comparable to the experimental value of 0.18 bar measured in this work. Three different AI potentials were used to verify the applicability of this approach. In contrast, when a classical force field (FF) was used to describe the interaction between CO2 and [P4444][CH3COO], the Henry’s law constant was computed to be 27 bar, significantly larger than the experimental value. The classical FF underestimates the CO2−[P4444][CH3COO] interaction compared with the AI calculations, which in turn leads to the smaller simulated CO2 solubility. Simulations further indicate that the CO2 interaction with the [CH3COO]− anion is much stronger than with the [P4444]+ cation. This result strongly suggests that the large CO2 solubility in [P4444][CH3COO] is due to the strong CO2−[CH3COO]− interaction.
1. INTRODUCTION Ionic liquids (ILs) are salts, which consist exclusively of cation and anion, with a melting point typically lower than 100 °C. ILs are attracting broad interests due to their remarkable properties, such as negligible volatility, high thermal stability, nonflammability, tunability, and solvation properties. ILs are of interest for a variety of gas-separation applications, but the most interest they received is in selective separations of CO2 from N2-containging streams in conventional power generation. This is due to stronger CO2 interaction with ILs than N2. For IL in solvent applications, CO2 solubility is the primary performance property of interest. Since the first report of high CO2 solubility in 1-butyl-3methylimidazolium acetate ([bmim] [CH3COO]) by Maginn et al.,1 subsequent experimental studies have shown high CO2 solubility in [emim][CH3COO] and [bmim][CH3COO] ILs.2−6 In addition, it has been shown that the [emim][CH3COO] IL exhibits a large CO2 permeability and a large permeability selectivity for CO2 over H2.7 Shiflett et al. have shown experimentally that CO2 forms complexes with [emim][CH 3 COO], [bmim][CH 3 COO], and [eeim][CH3COO].3,4,8 In their recent study,8 Shiflett et al. reported © 2014 American Chemical Society
the formation of imidazolium-2-carboxylates by way of a carbene intermediate instead of CO2 interaction with the basic [CH3COO]− anion proposed in their previous work.4 This CO2−carbene reaction mechanism has also been identified by Besnard et al.9−11 and Rogers and coworkers12 from their experimental studies with imidazolium acetate ILs. The [CH3COO]− anion is unique in that the CO2−carbene reaction occurs for the imidazolium acetate ILs; combinations of imidazolium cations with anions such as [Tf2N]−, [PF6]−, and [BF4]− do not show evidence of this type of reaction, and these ILs interact with CO2 through physisorption. One theoretical study performed by Cabaco and coworkers using ab initio (AI) gas-phase calculations concluded that CO2 is absorbed in imidazolium acetate ILs to form an imidazolium2-carboxylate through the CO2−carbene reaction mechanism.10 However, many other AI gas-phase calculations have shown that CO2 can interact strongly with the [CH3COO]− anion in the absence of the formation of a carbene.6,13 Received: March 10, 2014 Revised: May 15, 2014 Published: June 13, 2014 7383
dx.doi.org/10.1021/jp502425a | J. Phys. Chem. B 2014, 118, 7383−7394
The Journal of Physical Chemistry B
Article
In an effort to quantify the effect of the [CH3COO]− anion, we present here a method to calculate the CO2 Henry’s law constant for strong gas−solvent interaction systems by extending the test particle insertion method,20 which has been used to calculate CO2 Henry’s law constants in ILs15,21 through weak physical absorption. This goal is accomplished by explicitly including Henry’s law contributions from both linear and bent CO2 molecules using an AI-based potential. The method presented here is general and could be used for other systems. Additionally, we describe a procedure to obtain CO2− IL interactions from AI calculations. Developing a force field (FF) from AI calculations to describe CO2 interactions with ILs is an important but very challenging task. The first portion of the task is the generation of a database of CO2−IL interaction energies obtained from the AI calculations, which involves the selection of the AI theory and a procedure to account for the many-body effects in the condensed phase. The second portion of the task is the fitting of the database for later use in classical simulations. In our approach, we have constructed a “correction term” to be applied to a simple classical nonpolarizable force field. For simplicity, we do not consider the many-body effects here. In this work, the tetra-n-butylphosphonium acetate ([P4444][CH3COO]) IL system was chosen because the selection of [P4444]+ eliminates the CO2−carbene reaction path available in imidazolium containing ILs. The absence of this complication allows us to directly probe the CO2 interaction with the acetate anion both theoretically and experimentally. We first calculate CO2 interactions with both the [CH3COO]− and [P4444]+ ions in the gas phase using three different AI methods at many different CO2 positions relative to the anion and the cation. Once completed, the interaction results are used to obtain Boltzmann-averaged potentials, which can be used to calculate the Henry’s law constant. To the best of our knowledge, this is the first report of a computational technique to calculate the Henry’s law constant for gas absorption in solvents through strong interactions. Additionally, CO2 solubility in [P4444][CH3COO] was experimentally measured. The experimental Henry’s law constant was obtained from these measurements to validate the theoretical Henry’s law constant calculation for strong gas−IL interactions. In the following, we first describe the theoretical methods to compute the Henry’s law constant for physisorption and chemisorption systems, followed by the simulation and experimental details. The simulation results using a classical force field for physisorption and ab initio potential for both physisorption and chemisorption are presented along with the experimental measurements.
Recently, we have performed both AI gas-phase and AI condensed-phase calculations and found strong CO2 interactions with the [CH3COO]− anion.7 Both bent (Figure 1a)
Figure 1. Representative configurations for the [CH3COO]− anion interaction with bent CO2 (a) and linear CO2 (b) obtained from ab initio geometry optimization.7 Panel c shows the linear CO2−[P4444]+ interaction also obtained from ab initio geometry optimization. Linear CO2 interacts with ionic liquid through weak physisorption, and the ∠OCO angle is not significantly different from 180°.7 In contrast, bent CO2 interacts with ionic liquid through strong chemisorption, and the angles are significantly different from 180°.
and linear CO2 (Figure 1b) configurations have been identified to interact with [CH3COO]−.7,14 The linear CO2 interaction energy with [CH3COO]− is ∼10 kcal/mol, 6 kcal/mol stronger than with [Tf2N]− or [PF6]−. The interaction energy of bent CO2 with [CH3COO]− is 40 kcal/mol. All of the calculated strong CO 2 −[CH 3 COO] − interactions suggest that [CH3COO]− containing ILs should possess a large CO2 solubility. There are two types of methods to compute the Henry’s law constant, that is, the direct and indirect methods. In the direct method, one calculates the excess chemical potential, or solvation free energy of solute going from an ideal gas state to being dissolved in the solvent. The direct method includes the test particle insertion,15 expanded ensemble,15 transition matrix,16 and free-energy calculation17 methods. In the indirect method, one first calculates the absorption isotherm,18,19 then obtains the Henry’s law constant by linear fitting of the gas partial pressure versus the gas solubility on the isotherm at low gas solubilities. Among these methods, the test particle insertion one is the simplest. Additionally, all of these methods only account for physical absorption. In the case of CO2 Henry’s law constant calculation using the test particle insertion method, a fixed linear CO2 configuration from the gas phase (∠OCO = 180°) is randomly inserted in the solvent. The interaction energy between this inserted linear CO2 and solvent molecules is computed for subsequent Henry’s law constant calculation. Note that the linear CO2 configuration is used because it is the relevant configuration in the ideal gas phase. Consequently, these methods cannot be simply used to calculate Henry’s law constant for strong CO2−solvent interaction in chemical absorption, in which the CO 2 configuration is highly bent and the ∠OCO angle is significantly different from 180°.
2. THEORY 2.1. Henry’s Law Constant Calculation for Physisorbing Systems. In the NVT canonical ensemble, the chemical potential for one flexible CO2 molecule in NIL IL molecules at volume V and temperature T is given by
(
μCO ,liq = − k bT ln VqCO 2
2
∫ exp(−βUIL,CO ) 2
× r12 sin θ1 r22 sin θ2 drIL dr1 dθ1 dϕ1 dr2 dθ2 dϕ2 / 7384
( ∫ exp(−βU ) dr ) IL
IL
) (1)
dx.doi.org/10.1021/jp502425a | J. Phys. Chem. B 2014, 118, 7383−7394
The Journal of Physical Chemistry B
Article
where kb is the Boltzmann constant, β = (1/kbT), qCO2 denotes the partition function for the CO2 molecule corresponding to the kinetic energy for the three atoms of the CO2 molecule, and UIL,CO2 is the total potential energy for the solution system, which consists of NIL IL solvent molecules and one CO2 solute molecule. The UIL,CO2 includes the intra potential energies for CO2 itself, IL interactions among themselves (cation−cation, cation−anion, and anion−anion), and the CO2−IL interaction. The V and UIL are the volume and the potential energy for the neat IL solvent system. The rIL denotes the x, y, and z coordinates for the NIL molecules, and r1, θ1, ϕ1, r2, θ2, and ϕ2 correspond to the spherical polar coordinates for the solute CO2 molecule (Figure 2a). Note that the solution phase corresponds to a very low CO2 mole fraction of xCO2 = 1/(NIL + 1) ≃ 1/NIL for large NIL.
the Henry’s law constant can be obtained using the following equation H=
ρIL k bT exp( −βμex )
(4)
where NIL and exp( −βμex ) V ⎛ = ⎜ exp( −βUIL,CO2) × r12 sin θ1r22 sin θ2 drIL dr1 dθ1 ⎝
ρIL =
∫
⎞ ⎛ dϕ1 dr2 dθ2 dϕ2⎟ /⎜ ⎠ ⎝
⎞
∫ exp(−βUIL) drIL × Cintra⎠
⎟
Similarly, in the NPT ensemble the Henry’s law constant is also given by eq 4, where ρIL = (NIL/), is the average volume obtained from NPT molecular dynamics (MD) simulations, and the excess chemical potential μex is given by ⎛ exp(−βμex ) = ⎜ ⎝
∫ V exp(−βPV ) exp(−βUIL,CO ) × r12 2
⎞ sin θ1r22 sin θ2 dV drIL dr1 dθ1 dϕ1 dr2 dθ2 dϕ2⎟ ⎠ ⎛ /⎜ × ⎝
∫ exp(−βPV ) exp(−βUIL) dV
⎞ drIL × C intra ⎟ ⎠
Figure 2. (a) Schematic representation of CO2 using spherical polar coordinates. Points A and B are used to define θ1 and ϕ1. Note that the two points are set arbitrarily. The x, y, and z coordinates for the O1 atom (red); r1, θ1, and ϕ1 for the C atom (blue); and r2, θ2, and ϕ2 for the O2 (red) atom are used to determine the CO2 position and configuration. The distances r1 and r2 correspond to the C−O1 and C−O2 bond lengths, and θ2 corresponds to the ∠OCO angle. The angles are 0 ≤ θ1, θ2 ≤ π, 0 ≤ ϕ1, and ϕ2 ≤ 2π. (b) Schematic description of CO2 in a spherical shell. One O atom of the [CH3COO]− anion is located at the center of the sphere. The shell is bounded by r1 ≤ r ≤ r2, where r is the distance between the O atom (red) of the [CH3COO]− anion and the C atom (blue) of CO2. The orientation of the CO2 molecule is randomly generated and atom positions for [CH3COO]− are fixed.
=
2
k bTqCO C intra 2
PCO2,IG
(2)
where IG indicates the ideal gas. The Cintra is given by C intra =
∫ exp(−βUCO ,intra) × r12 sin θ1 r22 sin θ2 dr1 dθ1 2
dϕ1 dr2 dθ2 dϕ2
(5)
The μex in the NPT ensemble (eq 5) was calculated using the test particle insertion method as described: Step 1: NPT MD simulations were performed for the pure IL solvent using the classical FF, and snap shots were saved for analysis. Step 2: Hybrid Monte Carlo simulations (HMC)22,23 were performed for a system consisting of a large number of ideal gas CO2 molecules. Note that in the ideal gas there are no intermolecular interactions between different CO2 molecules. There are only intramolecular interactions between the atoms within the same CO2 molecule, such as bond and bend interactions. Step 3: One CO2 molecule was randomly picked from the ideal CO2 molecules obtained in Step 2. This CO2 molecule, with an internal coordinate of r1, r2, and θ2 (Figure 2a) determined in Step 2, was then rotated randomly for ϕ1 and ϕ2 both between −π and π and cos(θ1) between −1 and 1 to obtain the three Euler angles. Step 4: For each snapshot in Step 1, the CO2 molecule from Step 3 was inserted. The O1 atom of the CO2 molecule (Figure 2a) was randomly inserted in the simulation box, and the values for r1, r2, θ1, θ2, ϕ1, ϕ2 were obtained in Step 3. The interaction energy Uint between the inserted CO2 molecule and the solvent molecules was calculated using a classical FF. A value of V × exp(−βUint) was collected and finally averaged over the number of NPT snap shots and the number of CO2 insertions. Similarly, V was collected and averaged. The μex was computed using the second equation in eq 5. Note that in Step 2 a CO2 configuration in the range of r1 − r1 + dr1, r2 − r2 + dr2, and θ2 − dθ2 was generated with a probability proportional to exp(− βUCO2,intra) × r21 r22 sin θ2 dr1
The CO2 chemical potential in the ideal gas phase is given by μCO ,IG = −k bT ln
(3)
where UCO2,intra is the intra potential energy for a single CO2 molecule itself because there are no intermolecular interactions between CO2 molecules in the ideal gas phase. The UCO2,intra depends only on r1, r2, and θ2 and could be obtained from either a classical FF or from AI calculations. The Henry’s law constant is defined as H CO 2 = limPCO2,IG;xCO2→0(PCO2,IG/xCO2). By equating the CO2 chemical potentials in the solution (eq 1) and the ideal gas phase (eq 2), 7385
dx.doi.org/10.1021/jp502425a | J. Phys. Chem. B 2014, 118, 7383−7394
The Journal of Physical Chemistry B
Article
Henry’s law constant for bent CO2 to account for the fact that in the ideal gas phase the CO2 ghost molecule will adopt the linear-like configuration (∠OCO of 179° and CO bond length of 1.15 Å) for 2.9578 × 1018 times and adopts the bent-like configuration (∠OCO of 140° and CO bond length of 1.22 Å) for only one time. The large factor of w is due to the more positive energy of 27.48 kcal/mol for bent CO2 in the gas phase compared with the linear CO2, which was obtained from the ab initio gas-phase calculation at the MP2/cc-pVDZ level. At a different level of MP2/aug-cc-pVTZ, the w was calculated to be 8.8317 × 1019. To be consistent with the ab initio calculations to obtain the CO2−IL interactions at the MP2/cc-pVDZ level, the w value of 2.9578 × 1018 was used in this work.
dr2 dθ2. This was verified by comparing simulations and analytical calculations for CO2 using a classical FF (not shown here). The previously described procedure was used in this work to compute the Henry’s law constant for CO2 absorption in the IL when a classical FF was used to simulate the flexible CO2, the IL, and the CO2−IL interactions. 2.2. Henry’s Law Constant Calculation for Chemisorbing Systems. In principle, the above procedure could be used to directly compute the CO2 Henry’s law constant for the case of CO2 interaction with the IL through strong nonphysical interactions. However, the deficiency in this approach lies in the inability of a classical FF to accurately model both physisorption (weak binding) and chemisorption (strong binding) using one formalism and one set of parameters. An alternative method was developed and is described later. In our previous AI gas and AI condensed phase calculations,7 it is shown that CO2 could adopt both highly bent and roughly linear CO 2 configurations when interacting with the [CH3COO]− anion, which correspond to strong chemical and weak physical interactions, respectively. To simplify the Henry’s law constant calculation, thereby reducing the computational expense, we chose two specific CO2 configurations with θ(OCO) = 179°, r(CO) = 1.15 Å and θ (OCO) = 140°, r(CO) = 1.22 Å to represent the linear and bent CO2, respectively. Step 2 in the above approach was not used to generate CO2 configurations because a bent CO2 configuration corresponding to strong CO2−[CH3COO]− chemical interactions cannot be produced from the ideal gas-phase HMC simulation. Instead, a fixed rigid linear CO2 configuration (θ (OCO) = 179°, r(CO)=1.15 Å) was inserted in the simulation box to calculate the Henry’s law constant using the previously described Step 3 and Step 4. This will give the CO2 Henry’s law constant corresponding to physisorption. For the rigid bent CO2 (θ (OCO) = 140°, r(CO) = 1.22 Å), the Henry’s law constant was also computed in a similar way as the linear CO2, and this gives the Henry’s law constant for chemisorption. An ab-initio-based potential for CO2 interaction with IL corresponding to both physisorption and chemisorption was used in the calculation. Very importantly, in the case of bent CO2, the calculated Henry’s law constant must be scaled by a weighting factor of w when compared with the linear CO2. The weighting factor 1/w denotes the probability of bent CO2 relative to the linear CO2 in the ideal gas phase and is given by w=
3. SIMULATION DETAILS 3.1. Henry’s Law Constant for Physisorption Using Classical FF. The classical FF potential used to simulate IL, gas solutes, and IL−gas interactions is given by =(r ) =
k b(r − r0)2 +
bonds
+
∑
kθ(θ − θ0)2
angles
∑
kχ [1 + cos(n0χ − δ0)]
dihedrals
+
∑
k ψ (ψ − ψ0)2
impropers N−1
+
∑ i=1
⎧ ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ ⎫ qiqj ⎪ σij ⎥ ⎪ ⎢ σij ∑ ⎨4εij⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥ + ⎬ r rij ⎪ ⎝ rij ⎠ ⎦ j=i+1 ⎪ ⎭ ⎩ ⎣⎝ ij ⎠ N
(7)
where the symbols represent their conventional meanings.24 Standard Lorentz−Berthelot combining rules were used to calculate the mixed Lennard-Jones (LJ) interaction parameters. The LJ potential was switched from 10.5 to 12.0 Å. A Verlet neighbor list with a 13.5 Å radius was used. The intramolecular electrostatic and LJ interactions for atoms separated by exactly three consecutive bonds were scaled by 0.5 and were neglected for atoms separated by fewer than three consecutive bonds. The classical FF parameters for flexible CO2 molecule, the [CH3COO]− anion, and the [P4444]+ cation were obtained from the previous work.7,25,26 In the case of using a classical FF to describe the CO2, IL, and CO2−IL interactions, the simulation details of the test particle insertion method to calculate the Henry’s law constant are similar to those used in the work of Shah et al.15 Three sets of NPT MD simulations were performed starting from different initial configurations for pure [P4444][CH3COO] IL at 298 K and 1 bar by using the classical FF for the IL. The simulation box consisted of 100 ion pairs. Each simulation was typically run for 10 ns using the NAMD program,27 and 1000 equilibrated liquid structures corresponding to the last 2.5 ns were stored every 2.5 ps. To improve the sampling, we have used a similar “excluded volume” method15 that divides each stored structure into “occupied” and “unoccupied” regions. The cube length was roughly set to be ∼0.8 Å. Note that this cube length is slightly changed to give integer numbers of cubelets corresponding to the respective simulation box lengths. A cubelet was marked as occupied if the distance between the center of the cubelet and any IL atom was less than rev = s × (σi/2 + rprobe), where σi is the van der Waals parameter for the atom of the IL molecule, s is the scaling factor, and the radius
[exp( −βUCO2,intra) × r12r22 sin θ2]linear CO2 [exp( −βUCO2,intra) × r12r22 sin θ2]bent CO2
∑
(6)
The UCO2,intra was obtained from AI calculations. The value of w was calculated to be 2.9578 × 1018 at the MP2/cc-pVDZ level of theory. Note that eq 5 is the most general formula to calculate the Henry’s law constant for either physisorption or chemisorption. As shown in eq 5 and eq 3, the ghost CO2 molecule interacting with the solvent molecules in the test particle insertion method must be chosen from the CO2 ideal gas phase. This ghost CO2 molecule in the ideal gas phase with a set of r1, r2, and θ2 value (Figure 2a) exhibits a probability of exp(−βUCO2,intra) × r21 r22 sin θ2. In the case of using both linear and bent CO2 configurations as the ghost molecule, the probability of exp(−βUCO2,intra) × r21r22 sin θ2 for the ghost CO2 molecule must also be satisfied in the ideal gas phase. Hence, the w value was used to scale the 7386
dx.doi.org/10.1021/jp502425a | J. Phys. Chem. B 2014, 118, 7383−7394
The Journal of Physical Chemistry B
Article
3.2. Henry’s Law Constant for Physisorption and Chemisorption Using Ab Initio Potential. 3.2.1. Ab-InitioBased Boltzmann-Averaged CO2−IL Interaction Potential. The AI quantum single-point calculations in the gas phase were performed using the Gaussian 09 program.29 The interaction energy Uint,AI was calculated as the total energy of the dimer minus the energies of the two monomers at their respective geometries in the dimer system. The dimer was defined to be the CO2−ion cluster, where the ion could be either a cation or an anion. The monomers are CO2 and the ion. The calculations were carried out by randomly placing the C atom of a CO2 molecule in a spherical shell bounded between r1 and r2, as shown in Figure 2b for the CO2−[CH3COO]− dimer. The three Euler angles (θ1, ϕ1, and ϕ2) were randomly produced as in Step 3. Note that the two O atoms of the CO2 molecule are not necessarily located within the spherical shell. 20,000 acetate-CO2 configurations were generated for each spherical shell. The acetate anion and CO2 geometries were fixed; only the distance and orientation with respect to each other were varied. The frozen [CH3COO]− structure corresponded to the monomer configuration in the optimized bent CO 2 [CH3COO]− dimer in the gas phase. The frozen CO2 molecules were defined by values of θ(OCO) = 179°, r(CO) = 1.15 Å for linear CO2 and θ(OCO) = 140°, r(CO)=1.22 Å for bent CO2. Using this procedure, we generated CO2 configurations in spherical shells of 1.0 ≤ r ≤ 1.5 Å, 1.5 ≤ r ≤ 2.0 Å, 2.0 ≤ r ≤ 3.0 Å, 3.0 ≤ r ≤ 4.0 Å, 4.0 ≤ r ≤ 5.0 Å, 5.0 ≤ r ≤ 6.0 Å, 6.0 ≤ r ≤ 8.0 Å, and 8.0 ≤ r ≤ 10.0 Å. It has been shown that7 bent CO2 interacts strongly with the O atom of [CH3COO]− at a distance of ∼1.6 Å. Hence, two spherical shells with 1.0 ≤ r ≤ 1.5 Å and 1.5 ≤ r ≤ 2.0 Å between 1 and 2 Å were used to generate more data and better accuracy for Uint,AI. There were ∼1.6 × 105 such dimers for one set of frozen CO2 and [CH3COO]− configurations in the r range of 1 to 10 Å. Furthermore, to investigate the effect of the [CH3COO]− structure on Uint, another two sets of [CH3COO]− structures from NPT pure IL MD simulations were used. The total number of dimers for the CO2−[CH3COO]− system was 6 × 1.6 × 105, approximately 1 million. Because of the large number of AI single-point calculations, the MP2/cc-pVDZ level of theory was used instead of the more expensive MP2/aug-cc-pVTZ method, even though the latter is more accurate. To account for the difference between these two levels of theory, we performed single-point calculations for some of the dimers at the MP2/aug-cc-pVTZ level. Typically, the CO2 interaction with the [CH3COO]− anion at the MP2/ aug-cc-pVTZ level was found to be about 0.9 to 1.8 kcal/mol weaker than at the MP2/cc-pVDZ level. Despite this interaction energy difference, the major results obtained in this work were still found to be valid. Note that there were structures created in which the atom−atom distances between CO2 and IL were very short (0.06 to 1.0 Å). Because the interaction energies for such structures are large and unfavorable, such structures do not make contributions to exp(−βUint,AI). In each spherical shell, the value of exp(−βUint,AI) was
for the probe atom, rprobe, was set to be 1.0 Å. Note that a periodic boundary condition was used when calculating rev. When using the classical FF to describe the CO2−IL interactions, s was chosen to be 0.8. For each stored configuration, a list of cubelets was generated to record whether the cubelet is occupied (O) or unoccupied (U). The exclude volume fraction was calculated to be fev = (#U)/(#U + #O), where #U and #O denote the number of unoccupied and occupied cubelets, respectively. After the cubelet list for the configuration was generated, 2 million trial insertions were attempted to evaluate the excess chemical potential. Each trial insertion was performed in the following way. (1) The O1 atom of the CO2 molecule was randomly inserted in the box. If the O1 atom was in an occupied cubelet, this did not contribute to the excess chemical potential and a new CO2 molecule was inserted. If the O1 atom was in an unoccupied cubelet, the calculation continued to the next step. (2) A CO2 configuration was randomly chosen from the HMC simulation, as previously discussed in Steps 2 and 3 (Section 2.1). Using the O1 atom position in 1, the atomic positions for the whole CO2 molecule were generated. Note that the HMC simulation box for the ideal gas contained 500 CO2 molecules. About 50 000 HMC simulations were performed for equilibration. Every 100 steps during the 2million trial insertions, HMC was performed for the ideal simulation box. (3) The CO2 interaction with the IL molecules, Uint,FF, was calculated. The exp(−βUint,FF) was collected for excess chemical potential and subsequent Henry’s law constant calculations. Note that to speed up the calculations, we have used a modified Wolf’s method to calculate the electrostatic interaction between CO2 and all other IL molecules.25,28 The difference in CO2 electrostatic interaction with the IL was small between the modified Wolf’s method and the standard Ewald method. For 1000 different configurations, with each one consisting of one CO2 molecule and 100 ion pairs, it was found that the absolute average difference in CO2 electrostatic interaction with IL molecules between the modified Wolf’s method and the Ewald method was 0.26 kcal/mol, with the absolute difference ranging from 0 to 0.58 kcal/mol. After performing 2 million steps for the above one to three calculations and calculating the average of exp(−βUint,FF), it is very important to scale this value by a factor of fev to remove the biasing due to the use of the cubelet approach. When using eq 5 to calculate the Henry’s law constant, the O1 atom of the CO2 molecule should be inserted randomly in the whole simulation box rather than only in the unoccupied cubelets. If one inserts the O1 atom in the unoccupied cubelets to improve the sampling for n steps as used in the cubelet approach, this will implicitly mean that one has also attempted to insert the O1 atom in the occupied cubelets for n × #O/#U steps; even those insertions corresponding to the occupied cubelets do not contribute to the Henry’s law constant calculation due to the very repulsive interaction between CO2 and the IL molecules. The corresponding total number of O1 insertions would be n/ fev instead of n if one inserts the O1 atom randomly in the simulation box. Hence, the fev factor must be used to remove the bias in the cubelet approach. We have verified this scaling from our test calculations for a Lennard-Jones model system. The Henry’s law constant obtained from the simple insertion method without the use of cubelets was the same as the value obtained from the cubelet approach combined with the use of the scaling factor fev.
4
expcollected and an average was calculated as ((∑i=2×10 i= (−βUint,AI))/(2 × 104)). Note that the previously generated dimers correspond to the configurations on which exp(−βUint) are evaluated when using the test particle insertion method to compute the Henry’s law constant (eqs 4 and 5). Consequently, a Boltzmann-averaged interaction potential 7387
dx.doi.org/10.1021/jp502425a | J. Phys. Chem. B 2014, 118, 7383−7394
The Journal of Physical Chemistry B
Article
attractive CO2 interactions with [CH3COO]− at short distances. Hence, these regions that would be neglected when using a classical FF are explicitly included for consideration in the case of AI potential. Additionally, a fixed rigid CO2 configuration (linear or bent) was used in 2 instead of flexible CO2 configurations obtained from the HMC simulation. In this work, we only add the for the closest pairs of CO2−[CH3COO]− and CO2−[P4444]+. All other pair interactions for CO2−[CH3COO]− and CO2−[P4444]+ were calculated using the classical FF. The closest CO 2 − [CH3COO]− pair was defined as the [CH3COO]− anion, which exhibits the closest distance between the C of the CO2 molecule and the closest O of the acetate anion; while the closest CO2−[P4444]+ pair was defined as the [P4444]+ cation, which exhibits the closest distance between the phosphorus atom of the cation and the C atom of the CO2 molecule (Figure 1). If all pairs of CO2−IL interactions were corrected by adding , it was found that bent CO2 exhibits unreasonably strong interactions with all of the anions. This suggests that an effective pair potential for CO2−IL interaction by considering the many-body effect in the condensed phase should be used if all CO2−IL pair interactions are corrected. To investigate the effects on the CO 2 −[CH 3 COO] − interaction using different AI levels of theory, we also performed calculations at the BLYP/cc-pVDZ level using the Gaussian 09 program and M11 level using the Q-Chem.30 The M11 functional is a range-separated hybrid meta-GGA.31 This M11 DFT level of theory was previously shown to accurately reproduce gas-phase CCSD(T) interaction energies for a series of CO2−[CH3COO]− configurations.14 3.3. Calculation of Absorption Isotherm Using a Classical FF. For comparison, the absorption isotherm for CO2 in the IL was calculated using the classical FF. This was accomplished by applying the continuous fractional component (CFC) Monte Carlo (MC) method to compute the CO2 solubility in [P4444][CH3COO] at various pressures. The classical FF was used to describe the IL−IL, CO2−IL, and CO2−CO2 interactions. The simulation details of CFC MC simulations were very similar to those used in our previous work.25
bltz between CO2 and the acetate anion could be computed using the following equation
i = 2 × 104
exp( −β
)=
∑i = 1
exp( −βUint,AI)
2 × 104
(8)
In this way, was calculated for the CO2 interaction with the [CH3COO]− anion in all spherical shells. Note that the value of obtained from eq 8 is temperaturedependent. For example, was found to be 0.61 to 1.5 kcal/mol larger (more positive) at 373 K than at 298 K for the interaction of both linear and bent CO2 with [CH3COO]−. To be consistent, was computed at 298 K for all shells. Similarly, for all 2 × 104 dimers in each shell, the classical FF was used to obtain the CO2 interaction with the [CH3COO]− bltz anion, Uint,FF. A Boltzmann-averaged potential corresponding to the classical FF was obtained using a similar method as for . A correction for the Boltzmannaveraged potential was calculated using the following equation bltz bltz bltz = < Uint ,AI> − < Uint ,FF>
(9)
In the case of the CO2−[P4444] interactions using the MP2/ cc-pVDZ method, dimers were generated in a similar way as the CO2−[CH3COO]− system. The r is defined as the distance between the C atom of CO2 and the P atom of the cation. Note that 12 spherical shells were used, ranging from 1 to13 Å with a thickness of 1 Å for each shell. The [P4444]+ structure corresponds to that of the linear CO2−[P4444]+ optimized dimer. Because of a large number of 56 atoms for each CO2− [P4444]+ dimer, only 500 CO2 positions were randomly bltz generated in each shell. Values for , , and were obtained for both the linear and bent CO2 int,corr interaction with the [P4444]+ cation in each spherical shell. The linear and bent CO2 conformations interacting with [P4444]+ are the same as those used to determine the interaction with the [CH3COO]− anion. Note that even though the optimized structure for the CO2−[P4444]+ dimer in the AI gas-phase calculation predicts the linear CO2 structure (Figure 1c), the bent CO2 interactions with [P4444]+ were still calculated. This is because the bent CO2 interactions with both the [CH3COO]− anion and the [P4444]+ cation are needed when computing the Henry’s law constant. In this work, the many-body effect is not considered in developing the interaction potential. Some of our AI gas-phase calculations show that the CO2−[CH3COO]− interaction in a system containing two ILs does vary somewhat compared with the calculations for the system containing one [CH3COO]− anion. For linear-like CO2, the difference is 1 to 2 kcal/mol. For bent CO2, the difference is 0.3−6 kcal/mol. It is difficult to evaluate the many-body effects on Henry’s law constant calculations because the development of an ab-initio-based potential to include the many-body effect is very challenging. Despite this, the CO2−[CH3COO]− interaction energy obtained from AI calculations is more accurate than that from the classical FF. 3.2.2. Henry’s Law Constant Calculations. In the case of using the AI potential to describe the CO2−IL interaction, the same previously described procedure to calculate the Henry’s law constant as for the classical FF was used with the following exceptions. After calculating Uint,FF in 3, a value of was added to correct for the AI interaction. Note that smaller s values of 0.55, and 0.35 were used for linear and bent CO2 configurations, respectively, because the AI potential still gives +
4. EXPERIMENTAL DETAILS The 1H, 13C, and 31P NMR spectroscopic analyses were completed using a Bruker Avance III 400 MHz spectrometer. Fourier transform infrared (FTIR) spectra were collected on a Nicolet Spectrum 100 with an attenuated total reflectance (ATR) apparatus. Density measurements of the ILs were performed using a Micromeritics Accupyc II 1340 pycnometer. Water content was determined using a Metrohm 860 Karl Fisher (KF) Thermoprep titration unit, equipped with an 831 KF coulometer. Glacial acetic acid (2.12 g, 35.4 mmol) and tetrabutyl phosphonium hydroxide solution (40% aqueous solution, 25.1 g, 36.3 mmol) were stirred at room temperature in air for 5 h. The solvent was evaporated and the product was extracted into 60 mL of EtOAc, filtered, evaporated, and vacuum-dried to give a clear, viscous liquid. This liquid was vacuum-dried on a Kugelrohr at 110 °C for 5 days to give 11.5 g dried product (35.0 mmol, 99% yield). The dried IL contains no detectable water as it was subjected to Karl Fisher titration prior to CO2 absorption. C18H39O2P: FW = 318.47. density = 0.9703 g/mL at 23.8 °C. H2O content = 6160 ppm. Tdec = 325 °C. FTIR (ATR film): 2958, 2931, 2872, 1576, 1458, 1373, 1318, 1097, 7388
dx.doi.org/10.1021/jp502425a | J. Phys. Chem. B 2014, 118, 7383−7394
The Journal of Physical Chemistry B
Article
903, 719, 630 cm−1. 1H NMR (400 MHz, CDCl3): 2.31 (m, 8H, P-CH2), 1.88 (s, 3H acetate CH3), 1.47 (m, 16H, P-CH2), 0.92 ppm (m, 12H, P-CH3). 13C NMR (101 MHz, CDCl3): 176.8 (acetate CO2), 24.9 (acetate CH3), 23.6 (dd), 18.5 (d), 13.3 ppm. 31P NMR (162 MHz, CDCl3): 33.0. CO2 solubility measurements were performed using a PCTPro 2000 apparatus from Setaram. The PCT-Pro 2000 is a Sievert’s apparatus, which determines gas solubility in a sample by charging a sealed sample chamber of known volume with a known quantity of CO2. The CO2−charged sample chamber is isolated from the rest of the system, and the pressure drop in the chamber is measured. The drop in pressure due to CO2 absorption into the liquid is then measured, and the quantity of CO2 absorbed into the liquid is determined from an equation of state, in this case, the NIST Standard Reference Database.32,33 For all tests, 0.5 g of sample was loaded into the sample chamber. The sample was held under a dynamic vacuum for 4 h prior to starting a test. During this evacuation, the sample was heated to 30 °C and was stirred at 300 rpm. After evacuation, testing was initiated by dosing the sample chamber with a known amount of CO2, and the sample was allowed to equilibrate for at least 4 h. Dosing was repeated to obtain an isotherm over a pressure range of 0 to 10 bar. Throughout all tests, the sample was maintained at 30 °C and was stirred at 300 rpm.
pressures and a physisorbing interaction at higher pressures once the chemisorbing sites are occupied. The simulated CO2 solubility using a classical FF is small compared with the experimental data and does not exhibit chemisorption at low pressures. The largest CO2 solubility difference between the classical FF and the experimental data occurs at low pressures. These simulated CO2 solubility results indicate that the classical FF underestimates CO2 interactions with the IL, which is consistent with our previous AI calculations.7 As shown later, the classical FF underestimates the CO2−[P4444]+ and CO2− [CH3COO]− interactions compared with the ab initio calculations. This means that the simulated CO2 solubility values corresponding to the ab-initio-based CO2−IL potential will be larger than those solubilities using the classical FF. Finally, it is interesting to note that the experimental CO2 solubility values in [P4444][CH3COO] obtained in this work and [bmim][CH3COO] obtained by Shiflett3 are similar to each other. 5.2. CO2−IL AI Interaction Energy. 5.2.1. CO2− [CH3COO]− AI Interaction. MP2/cc-pVDZ Results. The interaction energy results for CO2 with one specific [CH3COO]− anion configuration at the MP2/cc-pVDZ level are shown in Table 1. Note that the Boltzmann-averaged interaction energies were rounded to two decimal digits, which typically corresponds to the simulated uncertainty. In the case of large simulation uncertainties, such as for bent CO2 interaction with [CH3COO]− at 3.0 ≤ r ≤ 4.0 Å, a large number of CO2 configurations (greater than 2 × 104) are needed to obtain better accuracy. Using the same [CH3COO]− anion as in Table 1, another 2 × 104 sets of different CO2 bent configurations were generated at 3.0 ≤ r ≤ 4.0 Å. The was found to be −29.48 kcal/mol, still reasonably close to the value of −29.96 kcal/mol (Table 1), even though the difference between these two values is in the first decimal place rather than in the second decimal place. For the linear CO2 interaction with [CH3COO]−, AI calculations show that is repulsive when CO2 is close to the [CH3COO]− anion at values of r ≤ 2.0 Å. The CO2 solubility contribution due to this region is small (small B value in Table 1). At large CO2−[CH3COO]− distances of r ≥ 6 Å, is weak, typically less than −0.5 kcal/mol. At 2.0 ≤ r ≤ 3.0 Å, the exhibits the strongest interaction of −7.20 kcal/mol, which corresponds to a CO2 solubility contribution of 60% (B(ri)/∑B(ri) = 60%) of the total linear CO2 solubility. In the case of bent CO2, the is still attractive at a small distance of 1.0 ≤ r ≤ 1.5 Å, but 14 kcal/mol weaker than the strongest interaction of −29.96 kcal/mol at 3.0 ≤ r ≤ 4.0 Å. At larger distances of r ≥ 6.0 Å, the bent CO2 interaction with [CH3COO]− is weak, approximately −1.0 kcal/mol. The CO2 solubility contribution exhibits the largest value at 3.0 ≤ r ≤ 4.0 Å, and significant CO2 solubility is also observed at values of 1.5 ≤ r ≤ 3.0 Å. Compared with linear CO2, the bent CO2 exhibits much stronger values (e.g., − 29.96 kcal/mol for bent CO2 compared to −7.2 kcal/mol for linear CO2). However, this large interaction energy difference does not necessarily imply that the CO2 solubility is dominated by the bent configuration because the large weighting factor w of 2.9875 × 1018 needs to be considered as well. When comparing the solubility of linear and bent CO2, it is interesting to note that bent CO2 shows high values of B (indicating significant solubility) over the range of 1.5 ≤ r ≤ 4.0 Å, while the linear CO2 exhibits high values of B over the range
5. RESULTS AND DISCUSSION 5.1. CO2 Absorption Isotherm. The experimental and simulated CO2 solubilities in [P4444][CH3COO] using a classical FF are shown in Figure 3. Experimentally, CO2 exhibits a significant solubility in [P4444][CH3COO] at low pressures, which starts to level off at high pressures. This behavior is typical of a chemisorbing interaction at low
Figure 3. Experimental (filled circles) and simulated (open squares and open triangles) CO2 solubility (mole fraction) in the [P4444][CH3COO] ionic liquid (IL) by using the classical force field (FF) for flexible CO2, IL, and CO2−IL interactions. The experimental CO2 solubility was measured at 303 K, and the simulated solubilities were obtained at 298 (squares) and 303 K (triangles) from the continuous fractional component Monte Carlo calculations. The experimental CO2 solubilities in [bmim][CH3COO] at 298 K obtained by Shiflett3 are also included for comparison (open diamonds). The simulated CO2 solubility is small compared with the experimental data because the classical FF underestimates the CO2−IL interactions. The experimental uncertainty was estimated to be 0.01 for CO2 solubility in [P4444][CH3COO] obtained in this work and 0.006 for CO2 in [bmim][CH3COO] obtained by Shiflett. The simulation error bars are typically less than the size of the symbols, and the experiment and simulation error bars are not shown for clarity. 7389
dx.doi.org/10.1021/jp502425a | J. Phys. Chem. B 2014, 118, 7383−7394
The Journal of Physical Chemistry B
Article
Table 1. Boltzmann-Averaged Interaction Potential (eq 8) at 298 K between CO2 (Linear and Bent Configurations) and One Specific [CH3COO]− Aniona AI (kcal/mol)
r (Å)
FF
A
1.0 1.5 2.0 3.0 4.0 5.0 6.0 8.0
≤ ≤ ≤ ≤ ≤ ≤ ≤ ≤
r r r r r r r r
≤ ≤ ≤ ≤ ≤ ≤ ≤ ≤
1.5 2.0 3.0 4.0 5.0 6.0 8.0 10.0
95.02 10.39 −7.20 −6.40 −5.36 −3.62 −0.40 −0.08
2.3 (22) × 10−70 2.4 (8) × 10−8 1.9 (3) × 105 4.9 (12) × 104 8.5 (7) × 103 4.5 (7) × 102 1.97 (3) 1.149 (3)
1.0 1.5 2.0 3.0 4.0 5.0 6.0 8.0
≤ ≤ ≤ ≤ ≤ ≤ ≤ ≤
r r r r r r r r
≤ ≤ ≤ ≤ ≤ ≤ ≤ ≤
1.5 2.0 3.0 4.0 5.0 6.0 8.0 10.0
−15.89 −29.86 −29.22 −29.96 −23.53 −8.63 −1.05 −0.23
1.5 (10) × 10−7 2.6 (12) × 103 880 (867) 3.1 (30) × 103 0.06 (4) 7.2 (36) × 10−13 2.0 (1) × 10−18 4.94 (2) × 10−19
B linear CO2 2.3 (22) × 10−69 5 (2) × 10−7 1.5 (3) × 107 7.7 (18) × 106 2.1 (2) × 106 1.7 (3) × 105 2.44 (4) × 103 2.349 (6) × 103 bent CO2 1.5 (10) × 10−6 5.0 (24) × 104 7.0 (69) × 104 4.8 (48) × 105 15 (9) 2.7 (14) × 10−10 2.5 (1) × 10−15 1.010 (4) × 10−15
(kcal/mol)
A
− 50.94 −2.88 −2.65 −2.17 −1.43 −0.27 −0.05
− 5 (4) × 10−38 128 (8) 87 (4) 39 (1) 11.1 (4) 1.59 (2) 1.087 (3)
− 17.27 −9.98 −8.98 −8.53 −5.97 −1.71 −0.46
− 7 (4) × 10−32 6.9 (7) × 10−12 1.3 (7) × −12 6 (1) × 10−13 8 (3) × 10−15 6.1 (6) × 10−18 7.37 (5) × 10−19
Results were obtained from ab initio (AI) calculations at the MP2/cc-pVDZ level. About 20 000 random CO2 positions relative to the [CH3COO]− anion were generated in each spherical shell. The r denotes the distance between the C atom of CO2 and the O atom of the [CH3COO]− anion 000 exp(−βUint,AI))/(20 000 × wCO2)) and B = A × 4/3π(r32 − r31) are also shown. The weighting factor w was (Figure 2b). Two values of A = ((∑i=20 i=1 set to be 1 for the linear CO2 (∠OCO = 179°, rCO = 1.15 Å) and 2.9578 × 1018 for the bent CO2 (∠OCO = 140°, rCO = 1.22 Å) molecule. The B value roughly corresponds to the CO2 solubility contribution in the respective shell. For comparison, the results obtained using the classical force field (FF) are also shown. Uncertainties in the last digit obtained from block average calculations are given in parentheses. The symbol “−” indicates that results were not computed. a
2.0 ≤ r ≤ 6.0 Å. This result suggests that CO2−[CH3COO]− interactions at shorter distances (r ≤ 2.0 Å) involve bent CO2 configurations, which could be interpreted as chemisorbing interactions. It has been shown experimentally that the CO2 solubility is small in ILs containing the [BF4]− anion.34−36 For example, the CO2 Henry’s law constant in [bmim][BF4] is 52.9 bar at 298 K,35 about 2300 times larger than the CO2 Henry’s law constant of 0.0222 bar3 in [bmim][CH3COO]. For comparison with CO2−[CH3COO]− interactions, AI calculations were − performed to compute between CO2 and the [BF4] anion in the same spherical shells as in Table 1. The interaction between linear CO2 and [BF4]− is typically 2 to 3.5 kcal/mol weaker than with the [CH3COO]− anion. In the case of bent CO2, the interaction with [BF4]− is significantly weaker than with the [CH3COO]− anion. For example, in the r range 1.5− − 4.0 Å, for the bent CO2−[BF4] interaction is 20 kcal/ mol weaker than the bent CO2−[CH3COO]− interaction. This weaker interaction is expected because the CO2−[BF4]− dimer optimization strongly favors the linear CO2 configuration. Accordingly, the bent CO2 contribution to CO2 solubility in [BF4]− should be negligibly small compared with the linear CO2 configuration. In both the linear and bent CO2 cases, CO2 interactions with the [CH3COO]− anion are stronger than with [BF4]−. This result is in agreement with the experimental observation that [bmim][CH 3 COO] has a higher CO 2 solubility than [bmim][BF4]. BLYP/cc-pVDZ and M11 Results. The CO2 interactions with the [CH3COO]− anion at the BLYP/cc-pVDZ and M11 levels of theory are shown in Figures 4 and 5. Compared with the MP2 calculations, the density functional theory of BLYP and M11 typically underestimates the linear CO2 interactions with
Figure 4. (a) Boltzmann-averaged interaction potential at 298 K between the linear CO2 (∠OCO = 179°, rCO = 1.15 Å) and the [CH3COO]− anion. The horizontal blue lines indicate the potential energy obtained from ab initio (AI) calculations at the MP2/cc-pVDZ level of theory in different spherical shells; the horizontal red lines correspond to the classical force field (FF) calculations in different spherical shells. The square and circle symbols correspond to AI calculations at the BLYP/cc-pVDZ and M11 levels of theory, respectively. The parameter r denotes the distance between the C atom of CO2 and the O atom of the [CH3COO]− anion (Figure 2b). (b) Boltzmann-averaged interaction potential at 298 K between linear CO2 and the [P4444]+ cation. The parameter r is the distance between the C atom of CO2 and the P atom of the [P4444]+ cation. The AI results correspond to MP2/cc-pVDZ calculations. In panels a and b, the dotted and dashed lines are added to guide the eye.
[CH3COO]− but overestimates the bent CO2 interactions with [CH3COO]−. In the case of the linear CO2 interaction with [CH3COO]−, the BLYP calculations give 0.1 to 1.1 kcal/mol 7390
dx.doi.org/10.1021/jp502425a | J. Phys. Chem. B 2014, 118, 7383−7394
The Journal of Physical Chemistry B
Article
stronger interaction energy than MP2 except at a large distance of 5.0 ≤ r ≤ 6.0 Å. The calculations performed at all three levels of theory show substantially stronger CO2 interaction energy when compared with the classical FF modeling. 5.2.2. CO2−[P4444]+ AI Interaction. For the CO2 interaction with [P4444]+ (Table 2), AI calculations show that both linear and bent CO2 interact with the cation repulsively within distances of r ≤ 4 Å. At larger distances of r ≥ 10.0 Å, the bent and linear CO2−[P4444]+ interactions are very weak, close to 0 kcal/mol. At intermediate distances, the linear CO2 interaction with the cation is typically
