Atomistic Simulation of Gas Uptake and Interface-Induced Disordering

Article Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

pubs.acs.org/JPCB

Atomistic Simulation of Gas Uptake and Interface-Induced Disordering in Solid Phases of an Organic Ionic Plastic Crystal Vinay S. Kandagal,* Fangfang Chen,* Jennifer M. Pringle, and Maria Forsyth Deakin University, Melbourne, Australia, Institute for Frontier Materials, VIC 3125, Australia

Downloaded via KAOHSIUNG MEDICAL UNIV on August 21, 2018 at 16:28:18 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Organic ionic plastic crystals (OIPCs) are a unique class of materials that exhibit a short-range disorder on the molecular level but are ordered at higher length scales. Recent experiments in our group have shown that the OIPC methyl(diethyl)isobutylphosphonium hexafluorophosphate ([P122i4][PF6]) exhibits a high ideal selectivity of 30 with respect to CO2 and N2 at 35 °C. Here, we employ classical molecular dynamics simulations for studying gas uptake in the OIPC [P122i4][PF6] at different temperatures. Both adsorption and absorption of the gases CO2, N2, O2, and CH4 were estimated using a gas/solid interface model. The observed trend in gas uptake was CO2 > CH4 > O2 > N2. The CO2 uptake was found to be dependent on both the OIPC structure and temperature. Owing to phase transitions and intermolecular motions, the solid gets disordered with increasing temperatures. The study finds that such disordering effects can also be caused by interface effects, which can enhance gas absorption. The results qualitatively confirmed that the OIPC can be effective in separating CO2 not only from N2, as per the previous experimental study, but also from O2 and CH4. However, strategies to improve the free volume in the material become important.

1. INTRODUCTION Human activities over the past years have perturbed the natural carbon cycle and led to a net accumulation of CO2 in the atmosphere.1 This is expected to cause a rise in the global temperature1 and subsequently affect our current climatic conditions. Therefore, to reduce its atmospheric concentration, CO2 needs to be actively removed from the atmospherea strategy termed as “negative emission”.2 Hence, the technologies that are aimed at capturing CO2 at point sources (such as power plants) are expected to play a crucial role. Aqueous amine-based gas separation methods are currently widely used for large-scale CO2 separation in industries.3−5 A considerable amount of energy, however, is required to regenerate the amine-based solvent, and hence, this increases the cost of the entire process.3 Another problem in this process is the loss of solvent in the effluent gas.6,7 In this regard, a more energy efficient mechanism based on physical separation is attractive. Membrane-based gas separation can be such an alternative method and has been used in industrial gas separation applications for a long time.5,8,9 Polymers are the most widely used membrane materials for gas separation applications.9 Polymers above their glasstransition temperature (Tg) can exhibit high gas permeabilities because of large fractional free volume.10 However, for industrial applications, glassy polymers, which are below Tg, are used because of their better mechanical strength.10 Gas sorption in polymers is mainly driven by free volume; however, their chemistry can be tailored (e.g., attaching specific functional groups) to enhance selectivity based on gas solubility.11 © XXXX American Chemical Society

Recently, ionic liquids (ILs) have generated considerable interest as potential materials for CO2 capture.12−15 ILs are a relatively novel class of materials that have attracted considerable attention, mainly because of their nonvolatile nature and chemical tuneability.14,15 Both experiments and simulations have been carried out to understand the CO2 solubility in various ILs, which were prepared by a combination of different cations and anions.13,16−18 This knowledge can be relevant to the study of organic ionic plastic crystals (OIPCs) for gas separation as both classes of materials are made up entirely of cations and anions, with OIPCs as the solid-state counterparts of ILs. CO2 solubility in various ILs was found to be strongly influenced by the type of the anion19 and was further attributed to the strength of gas−anion interactions.20 However, Kazarian et al.21 noted that in certain CO2−IL systems, stronger anion− CO2 interactions corresponded to lower gas solubility. On the basis of this trend, the authors proposed that the anion−gas interactions solely did not dictate CO2 solubility in those ILs, and free volume also played an important role. Subsequent studies in the literature have indicated toward a greater role of free volume in influencing gas solubility in ILs.22−25 The ILs in membrane applications are used in combination with a solid polymer membrane, which provides mechanical stability, to form supported IL membranes (SILMs). Scovazzo26 Received: June 6, 2018 Revised: July 27, 2018

A

DOI: 10.1021/acs.jpcb.8b05444 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B noted that the reported ideal selectivity values for CO2/N2 in the literature for SILMs were in the range 15−57. In spite of their superior performance, SILMs can still suffer from loss of solvent at high and moderate pressures.27 OIPCs can help overcome this problem by their inherent solid state at room temperature. As they have similar chemical compositions to ILs, OIPCs can also be expected to exhibit good CO2/N2 selectivity. Recently, a novel membrane composed of polyvinylidene fluoride (PVDF) and an OIPC, methyl(diethyl)isobutyl phosphonium hexafluorophosphate ([P122i4][PF6]), was found to exhibit an ideal selectivity of ∼30 at 35 °C for CO2/N2.28 The OIPC [P122i4][PF6] exhibits three solid−solid-phase transitions, at 303, 353, and 393 K, before it starts to melt at around 418 K.29 Multiple solid−solid-phase transitions are commonly seen in OIPCs because of relatively weak cation− anion interactions, which make the temperature-induced molecular motion easier. A previous molecular dynamics (MD) study on [P122i4][PF6] identified different extents of rotational and translational motions of the ions at different temperatures, which can be linked to multiple solid−solid-phase transitions.30 In the current study, we use MD simulations to investigate the potential gas separation capability of this OIPC. The study of OIPCs for gas absorption is still in its infancy, and the fundamental understanding of the gas sorption mechanism and the impact of multiple solid−solid-phase transitions is still lacking. Therefore, in the current study, the adsorption and solubility of CO2, N2, O2, and CH4 in [P122i4][PF6] was investigated across different solid-state phases. Following the experimental study by McDonald et al.,28 the initial MD calculations31 on CO2 and N2 absorption in the OIPC demonstrated that the surface effects are crucial in modeling the OIPC solid. In that study, the solid OIPC model with vacuum interface became highly disordered at 325 K, which in turn had a positive effect on CO2 uptake. Such a structural disordering change has also been observed experimentally in experiments involving OIPC/PVDF composites.32 Hence, in the current study, we extend the MD simulation box to study a longer solid model, in the direction perpendicular to the interface, in order to further investigate the effect of the vacuum/ OIPC interface and the chemical nature of the permeating gas on OIPC disordering and gas solubility in [P122i4][PF6].

Figure 1. Schematic of the simulation model employed in the current study. The simulation box with relevant dimensions is shown along with a magnified view of the molecular structures of the [PF6] anion and the [P122i4] cation.

The MD simulation procedure used is similar to that used in our earlier work.31 Before introducing the vacuum and gas interfaces, the OIPC was initially simulated in bulk at four temperatures of 275, 310, 325, and 365 K, corresponding to its three different phases of IV (at 275 K), III (310 and 325 K), and II (365 K), as seen from the experimental analysis.29 The solid system was heated slowly from 123 K to these selected temperatures by a temperature step of 50 K. At each step, the temperature was ramped up in a time interval of 20 ps in the NVE ensemble followed by a slight relaxation for about 40 ps in the NPT ensemble. At the four specified temperatures, a longer NPT equilibration calculation was carried out for 4−10 ns. The vacuum region was then introduced, and the whole system was equilibrated again for different intervals of time at different temperatures until no further structural changes were observed in the solid. Finally, the CO2, N2, O2, and CH4 gases were introduced in the vacuum regions and allowed to permeate the solid. Both CO2 and N2 were introduced at an initial pressure of 26 and 40 atm, whereas O2 and CH4 were introduced at 40 atm. Such high pressures were used to ensure quicker gas permeation as the OIPC solid is considerably longer (12−13 nm) in the direction perpendicular to the gas/solid interface. The gas/solid system was then simulated in an NVT ensemble for different intervals of time depending on the temperature, the longest being 140 ns for the case of 310 K. While CO2 absorption was studied at all the four temperatures, N2 was studied at 325 and 365 K, and O2 and CH4 were simulated at 325 K. CO2, N2, and O2 were modeled as rigid molecules and treated according to the method by Kamberaj et al.33 The potential parameters for both CO2 and N2 were obtained from Potoff and Siepmann.34 These parameters have been used in the literature for gas absorption simulations in ILs.18,35,36 The parameters for O2 and CH4 were borrowed from Vujić and Lyubartsev37 and Sun et al.,38 respectively. All the gas molecules except methane were modeled as rigid molecules. The TraPPe force field for CO2 used in the current study was originally formulated as a

2. METHODOLOGY The schematic of the gas/solid interface molecular model employed in the current study is shown in Figure 1. The supercell of the OIPC solid consists of 2 × 2 × 8 unit cells, with each unit cell consisting of eight pairs of cations and anions. Its crystalline structure was obtained from a previous X-ray crystallography measurement carried out at 123 K.29 The vacuum regions were placed on both sides of the OIPC solid. Gas molecules were then introduced in these regions after initial equilibration of the solid with the vacuum. The solid was 12−13 nm longer along the Z direction (depending on the temperature), whereas the total length of the vacuum/gas regions was at least 15 nm. This was chosen to ensure the minimum influence of the periodic images of the OIPC atoms in the Z direction. The schematic in Figure 1 also shows the molecular structure of the [P122i4] cation and the [PF6] anion. The structural analysis presented in the Results and Discussion section mainly involves phosphorous atoms on both the anion and cation, denoted as P1 and P2, respectively. B


Article

The Journal of Physical Chemistry B

Figure 2. Number density profiles of the cation and anion phosphorous atoms (labeled as P2 and P1, respectively), at different temperatures, with an OIPC/vacuum interface. Only the right half of the solid is shown as the other half is similar. The horizontal dashed line represents the peak height in the bulk solid (i.e., without the interface).

rigid gas model. Both flexible and rigid model versions of this force field for CO2 have been used in the literature for the simulation of CO2/IL systems.18,35,39 From the literature, it can be seen that this rigid CO2 model was able to simulate the thermodynamic properties of CO2 in IL and polymer systems with good accuracy. For example, Huang et al.40 found a close match between simulation and experiment (with ∼2% deviation) for the partial molar volume of CO2 in the IL 1-nbutyl-3-methylimidazolium hexafluorophosphate ([bmim][PF6]) using the rigid CO2 model. Using the same, Hong and Panagiotopoulos41 found a good agreement between simulated and experimentally determined Henry’s constants for CO2 in polyethylene oligomers over a range of temperatures. Hence, the rigid CO2 model was employed in the current study. Similarly, rigid models for N2 and O2 were used, assuming that the effect of the bond vibrational motion is negligible in determining their solubility. The N2 model used in the current study was able to reproduce the right trend in the temperature dependence of solubility as per simulations carried out by Liu et al.42 in certain IL systems. The O2 parameters by Vujić and Lyubartsev37 were refined to match the experimental adsorption studies of the gas with zeolite systems. The rigid model treatment of the gas molecules is considered sufficient in obtaining a qualitative trend in their absorption in the OIPCs, which is the aim of the current study. One may incorporate the internal degrees of freedom (bond stretching and angle bending) for a more rigorous study involving comparison with the experimental data. The force field parameters for the OIPC solid used in the current study were formulated and tested previously by Chen et al.30,43 The program large-scale atomistic and molecular massively parallel simulator or LAMMPS (version 7 Dec 2015) was used for all MD calculations.44 A cutoff of 1.2 nm was used for both van der Waals and electrostatic energies in real space. The longrange electrostatic energies were treated using a particle− particle particle−mesh technique as implemented in LAMMPS. A time step of 1 fs was used for the initial 30 ns in all calculations, after which a time step of 2 fs was used for the further 60−110 ns calculations at 325 and 310 K. The SHAKE algorithm was used

to constrain all hydrogen atoms bonded to carbon while using a higher time step. To track the structural changes in the solid brought about because of interface effects, the atomic number density was calculated as a function of the Z coordinate using the standard method of histogram analysis using a bin width of 0.1 Å. Two-dimensional (2D) density distribution functions were similarly calculated, using codes developed in-house, to analyze the gas/ion structural information. Gas uptake involves both surface adsorption and absorption into the solid bulk. In the current work, absorption, which is a bulk property, is represented by defining an average density of solid the absorbed gas (e.g., CO2) in the solid (denoted as ρCO ̅ ). 2

gas Similarly, ρCO ̅ denotes the CO2 density in the bulk gas phase. 2

The gas adsorption at the gas/solid interface was calculated as per eq 1 as surface excess,45 and the total gas uptake, denoted as uptake , was calculated using eq 2. The interface boundary is NCO 2 positioned at the origin in the equations. They are described in detail in the Supporting Information. surf NCO = 2

0

∫−∞ (ρCOgas

2

gas gas − ρCO ̅ 2 )dV +

dV solid

∫0

∞

solid solid (ρCO − ρCO ̅ 2) 2

(1)

uptake surf solid solid NCO = NCO + ρCO ̅ V 2 2 2

(2)

surf In the above equations, NCO is the number of adsorbed gas 2 solid solid gives the total number of absorbed molecules, and ρCO ̅ V 2

uptake gas molecules. The total gas uptake, NCO , is therefore the sum 2 solid gas of these two terms. ρCO ̅ 2 represent the average CO2 ̅ 2 and ρCO density in the gas phase and inside the OIPCs, respectively. Vgas and Vsolid are the volumes occupied by the gas and the solid, respectively, inside the simulation box. The abovementioned method of analysis was similarly applied to all other gases.

C


Article


Figure 3. Snapshots from the MD simulation of the OIPC/vacuum interface along the XZ and YZ planes and at two different temperatures. The unfolded atomic coordinates over a simulation run time of 50 ns for 310 K and 25 ns for 325 K are shown.

3. RESULTS AND DISCUSSIONS To understand the effect of having an interface on the structure and dynamics of the OIPC, the initial discussion will focus on how the vacuum interface perturbs the ordering, in particular near the interface of the OIPC, and how far any effects permeate into the OIPC structure. This will have implications for the bulk absorption of CO2 and other gases considered in this work. The interactions between the OIPC and CO2 are mainly analyzed in different phases (different temperatures) to determine the effect of gas chemistry on the solubility and also which species of the OIPC (anion or cation) contribute most to the solubility of the gases. 3.1. Effect of the Vacuum Interface on the OIPC Structure. To track the structural changes, the local density profiles, as shown in Figure 2, were calculated for the phosphorous atoms on both the cation (P2) and the anion (P1), along the Z axis. As the density profile is symmetric about the center of the simulation box, only the right half of the profile is presented here. Large changes in these profiles are evident closer to the solid/vacuum interface, shown as the shorter and broader peaks compared to the taller peaks in the central portion of the solid. This indicates the presence of a greater translational disordering of the atoms near the interface. The horizontal dashed lines marked in Figure 2 represent the normal peak heights of the density profiles from the bulk phase simulation, that is, in the absence of a vacuum interface, obtained at the same temperature. By comparing it with the corresponding peak heights in the interface model, one can get an idea of the extent to which the interface affects the overall OIPC structure at different temperatures. At 275 and 310 K, the solid/vacuum system was equilibrated for 15 and 40 ns, respectively. At both temperatures, the vacuum interface mainly affects the OIPC structure closer to the interface, in a region with the Z length beyond about 30 Å from the center of the solid as shown in Figure 2a,b, as suggested by the decrease in the density peaks. The OIPC structure in the center portion of the solid, however, appears to be well-maintained, as the heights of the density profile match the dashed lines. Therefore, a long-range order was maintained in the solid at these temperatures. When the temperature was increased to 325 and 365 K, the interface affected the entire OIPC structure. The simulation times were 25 and 12 ns at 325 and 365 K, respectively. All the density peaks at these two temperatures were lower than the marked dashed lines obtained from the bulk-phase simulation. The density profile is almost flattened at 365 K, which is more representative of that of an IL,46 indicating that the OIPC structure became liquid-like at this temperature. It should be noted that

experimentally, the OIPC is in a solid phase II at this temperature, wherein the lattice structure exhibits a certain degree of order.29,30 Hence, in the simulation, the interface effects were large enough to increase the ionic motions and disorder the entire OIPC structure. Such interface effects on ionic motions should be decreased in the bulk region if a longer simulation box in the Z-direction was used, but this approach prohibitively increases the computational cost. We note that this has implications even for real experimental systems, where an interface between an OIPC and another surface may lead to different degrees of disordering (or indeed ordering). This is a topic also under investigation by us and will be presented in a future paper. Figure 3 shows the snapshots of the OIPC viewed in the XZ and YZ planes. The snapshots show unfolded coordinates (without periodic boundary conditions) of the molecules over a simulation time of 40 and 25 ns, respectively. The OIPC structure is relatively more disordered at 325 K. It is interesting to see from the literature that ILs showed an opposite behavior, wherein their structures were slightly more organized at the IL/ vacuum interface than in the bulk.47 Furthermore, the figure also shows that at 325 K, there was a relatively higher disordering along the Y direction compared to the X direction. This may be due to more free space in the Y direction surrounding an ion in the lattice. 3.2. Gas Adsorption and Absorption. After the OIPC solid/vacuum system was sufficiently equilibrated, gas molecules were introduced into the vacuum regions of the simulation box. The CO2 absorption was studied at all the aforementioned temperatures, whereas the N2 absorption was only studied at the higher temperatures of 325 and 365 K, involving a more disordered OIPC. In the case of ILs, the gases O2 and CH4 have been found to exhibit low solubility.15,19 Hence, expecting a similar behavior in our current study, these gases were studied at a higher pressure of 40 atm (for quicker equilibration) and at 325 K. The number of adsorbed and absorbed gas molecules was calculated using the equations described in the Methodology section. Both these quantities are expressed in terms of surface and bulk concentrations in Table 1. The surface concentration (or surface excess), Γ, is obtained as the number of adsorbed gas molecules (as per eq 1) per unit interfacial area. At 275 K, the gas molecules were mostly retained on the surface (shown in the Supporting Information as Figure S2) with none entering the solid bulk. The highest CO2 adsorption was also seen at this temperature. With an increase in the temperature from 275 to 365 K, the gas adsorption gradually decreased for CO2. The low bulk solubility at 275 K can be attributed to the fact that the D


Article

The Journal of Physical Chemistry B Table 1. (a) Surface and Bulk Concentrations, Denoted as Γ and ρ̅solid gas , Respectively, of CO2 and N2 in the OIPC Calculated at Different Temperatures and 26 atm Initial Gas Pressure; for N2, Calculations Were Carried Out at Only Two Temperatures and at 26 atm; (b) Same Calculations Carried Out by Including O2 and CH4 Gases, at 325 K and 40 atm Initial Gas Pressure, and the Corresponding Uncertaintiesa

the lack of sufficient free volume. This behavior is also consistent with our previous calculations,31 wherein a shorter solid model, with Lz = 6.3 nm, was used. In these simulations, carried out at 275 K, the gas was found to be mainly adsorbed onto the solid surface. Upon increasing the temperature to 310 K, a considerable increase in CO2 absorption was observed. The concentration of CO2 inside the solid reached as high as 345 mol/m3. However, it dropped to 224 ± 5 mol/m3 when the temperature was further increased to 325 K, and finally remained similar at 222 ± 2 mol/ m3 at 365 K. As discussed earlier, the jump in gas absorption is likely attributed to an increase in the free volume at 310 K. The temperature effects on CO2 solubility in ILs have been investigated in the literature. It has been found to decrease with an increase in the temperature in a number of ILs.15,19,48 Similar trends for gas absorption are also seen in polymers such as polystyrene, polypropylene, polyethylene, and poly(vinyl acetate).49,50 One might therefore expect the same behavior in the current OIPC. However, while the surface concentration decreased with an increase in temperature, the gas solubility increased considerably from 275 to 310 K, and then decreased, as expected, at higher temperatures. The jump in solubility from 275 to 310 K is attributed to the insufficient free volume in the solid at low temperatures, whereas, at 310 K, the subsequent expansion of the solid accompanied by disordering led to a higher CO2 solubility. Therefore, in the case of higher simulated

(a) N2

CO2 T (K)

Γ (μmol/m2)

275 310 325 365

3.6 ± 0.04 1.6 ± 0.07 1.4 ± 0.06 0.8 ± 0.2

3 ρ̅solid gas (mol/m ) 0 345 ± 10 224 ± 5 222 ± 2 (b)

gas (325 K)

Γ (μmol/m2)

CO2 N2 CH4 O2

5.99 ± 0.04 3.8 ± 0.06 2.5 ± 0.1 1.7 ± 0.1

Γ (μmol/m2)

0.4 ± 0.09 0.5 ± 0.08

3 ρ̅solid gas (mol/m ) NA NA NA NA

3 ρ̅solid gas (mol/m ) 287 ± 8 NA 105 ± 9 62 ± 4

a

NA refers to values not accessed in the simulations.

OIPC structure is still highly ordered and dense at this temperature, and thus making absorption difficult because of

Figure 4. Time evolution of density profiles at 310 K, showing the structural changes resulting from CO2 absorption. The red and blue peaks represent the P1 and P2 atomic densities, respectively, whereas the green plot represents the densities of CO2 carbon. E


Article


Figure 5. 2D contour plots of the atomic densities of phosphorous in the cation and anions (shown as P2 and P1 respectively), and in the XY plane. The plots were calculated in a central section of the OIPC at the initial 5 ns (a) and between 90 and 95 ns (b), at 310 K. The color bar represents densities in units Å−3.

Figure 6. Radial distribution plots between the CO2 atoms and different atoms in the cations and anions at three temperatures. The RDFs for X−CCO2 pairs are shifted upward by 2 units along the Y axis for clarity.

temperatures (310, 325, and 365 K), when the solid model had sufficient free volume to accommodate the CO2 molecules, its solubility was dictated by thermodynamic temperature effects observed in other materials as discussed above. The gas absorption simulations were carried out at a fixed volume (in the NVT ensemble); however, the solid was allowed to expand along the Z direction, into the vacuum regions. Interestingly, it was found that the solid volume at 310 K increased significantly with CO2 absorption and was ca. 3% higher than the solid volume at 325 K after CO2 absorption. This could be the reason for a higher bulk solubility of CO2 at 310 K because of a large volume expansion accompanied by disordering upon gas absorption. The volume changes at 325 and 365 K, however, were negligible during absorption. The differences observed in the different phases may be a result of the fact that the solid models at these higher temperatures were already partially disordered and thus had the free volume required to accommodate the additional gas molecules. In the case of N2, no gas molecules were found inside the solid region beyond a depth of 6 Å from the interface boundaries. Hence, the total gas uptake by the solid was attributed to surface excess. There was a negligible change in the surface excess of N2

with temperature. Furthermore, it exhibited negligible bulk solubility at high gas pressure. The CH4 and O2 gases showed lower absorption than CO2 in the bulk at a high pressure of 40 atm, at 325 K. The simulations, therefore, showed that among the gases considered, CO2 exhibited the highest uptake by the OIPC, followed by CH4, O2, and N2. Hence, the results suggest that the OIPC can be effective in separating CO2 from these other gases. This trend in gas uptake was also previously observed for imidazolium-based ILs.15,19,48 The considerable difference between CO2 and N2 uptake qualitatively agrees with the previous experimental finding of high selectivity in the material.28 The solid showed a higher CO2 absorption at 310 K than at 325 K, attributed to considerable disordering and volume expansion that accompany the absorption at the lower temperature. Total simulation times at both temperatures were 146 and 98 ns for 310 and 325 K, respectively. The simulations were repeated with a different initial structure at 310 K to double check this phenomenon, and the result was found to be consistent. Figure 4 shows the structural changes in the OIPC as a result of gas absorption, in terms of atomic density peaks obtained at F


Article


Figure 7. (a,b) 2D density distribution plots for CO2−P1 (anion) and CO2−P2 (cation) interactions, (c) snapshot from MD simulation showing the cation and anion molecules within a radius of 5 Å from a CO2 molecule, and (d) coordination number of the P atoms around a CO2 molecule. The color bar represents densities in units Å−3. All the plots are obtained at 325 K.

one parallel to and the other perpendicular to the CO bonds. This gives a consistent view of the environment of the molecule along these two vectors, denoted as h and r in the plots. Atomic densities in the 2D plots were evaluated using eq 3.52 Therefore, the term ⟨N(h,r)⟩ in the equation represents the ensembleaveraged number of the neighboring atoms at a point (h,r) from the CO2 carbon atom. A spatial resolution of 0.2 Å, denoted by Δh and Δr in the equation, and a common cutoff of 8 Å along both the r and h axes was used. As the RDFs were similar across different temperatures, the 2D plots are shown only for the case of 325 K. Furthermore, the plots include densities of ions averaged over all CO2 molecules in the system, that is, both at the surfaces and inside the OIPC solid.

six different periods of simulation time. The peak heights of the P atoms in the cation and anion of the solid reduced gradually toward the end of the simulation, with greater reduction at the interfaces. This reduction is accompanied by ∼3% volume expansion in the solid, which was noted by tracking the shift in the peak positions at the interfaces marked by two dotted lines. 2D density contour maps were generated for the phosphorous atomic densities for the first 5 ns and between 90 and 95 ns in the simulation and are shown in Figure 5. The plots were obtained over the central region of the OIPC with a volume equal to that of a unit cell. A unit cell of the OIPC contains eight pairs of anions and cations, as shown in the plots. The plots demonstrate the disordering effect with time, as ionic densities that were initially localized became more dispersed. 3.3. Determination of the Nature of the CO2−Ionic Interactions from the Analysis of RDFs and 2D Density Plots. Figure 6 shows the radial distribution function (RDF) plots of CO2 interacting with various atoms in the cations and anions at the three temperatures 310, 325, and 365 K. For all the pairs, the positions of the first peaks on the X axis coincide at all the three temperatures, suggesting similar gas−ionic interactions across all the phases. The first CO2−P1 (anion) peak appears at 4.1 Å and is similar to that predicted by simulations for ILs consisting of the same anion.18,51 The first minimum in the RDFs is taken as the first solvation shell. Hence, the first solvation shell of CO2 is predominantly made up of [PF6]− anions and the alkyl group carbons in the cations (as suggested mainly by the C−OCO2 RDF plots). Our previous work showed that CO2 interacted closely with fluorine atoms in the anions and hydrogen atoms in the cations.31 This is also shown in Figure 6db as the H−OCO2 and F−CCO2 pairs peak before other pairs. 2D spatial distribution plots were also generated as shown in Figure 7, to gain insights into the environment surrounding an absorbed CO2 molecule. The plots present the probability density of finding cations and anions within a cylindrical region around the gas molecules. The distance vector between the CO2 carbon and a nearby atom was resolved into two components,

ρ (h , r ) =

⟨N (h , r )⟩ 2πr ΔhΔr

(3)

Figure 7a shows the density distributions for P1 atoms (in the anion) within an 8 Å radius from the CO2 carbon. The color bar represents the density values in the distribution. The plots indicate the presence of multiple anion locations around a CO2 molecule, within the considered cutoff. The anions are located in a direction perpendicular to the CO bonds and are closer to the CO2 carbon. Similarly, Figure 7b shows the density of the cation phosphorous atoms (P2), which are mostly located closer to the oxygen atoms. This is due to the Coulombic attraction between the negatively charged oxygen and the positively charged P2 atoms. The information in both the plots is summarized in the snapshot shown in Figure 7c. For clarity, a shorter cutoff of 5 Å was used in Figure 7c. The 2D plot of the density distributions is consistent with our previous results of ab initio calculations and the classical MD simulations.31 The ab initio calculations showed that the optimum configurations consisted of a CO2 oxygen atom pointing toward the phosphorous atom of the cation, and the CO2 carbon being closer to the anion phosphorous atom. This was further confirmed by three-dimensional density distribution plots obtained from classical MD simulations in the same work. G


Article


Figure 8. Radial distribution plots between the phosphorous atoms of the cations and anions in the bulk ordered (a) and disordered structure (b), both at 325 K.

Figure 9. 2D distribution plots for the P atoms in anions (a) and cations (b) around a N2 molecule, obtained at 325 K and 26 atm initial gas pressure.

DAA′, and DBB′ coinciding with a minimum for the P2−P2 distribution. Hence, this suggests that there is a higher fraction of the required interatomic separation of molecules to accommodate CO2 in the disordered solid than in the ordered structure. In other words, the completely ordered OIPC structure does not offer a sufficient free volume for CO2 absorption. This requires that a certain degree of disordering be introduced in the OIPC to enhance absorption. 2D distribution plots similar to those in Figure 7 are obtained for the case of N2 and the P atoms as shown in Figure 9. It can be seen from the figure that the intermolecular distances are similar for both the gases when interacting with the P atoms of the OIPC ions. Similar plots generated for all the cation carbon atoms around the gas molecules are shown in the Supporting Information as Figure S3, which showed that the carbons are slightly closer to N2 than CO2. However, in all these cases, the atomic densities around N2 are considerably lower than those found around the CO2 molecules in Figure 7a,b, as denoted by the scale of the color bars. This suggests that CO2 exhibits stronger interactions with the OIPC atoms than N2 under similar conditions. This difference in interaction strengths has led to lower adsorption for N2 than CO2 as noted in Table 1a. Furthermore, our previous ab initio-based calculations31 have also shown that a N2 molecule interacting with a single [PF6] or [P122i4] molecule in a vacuum exhibited lower binding energies compared to the corresponding values obtained for CO2. Hence, the stronger interaction of the OIPC ions with the more polar CO2 gas rather than free volume seems to be a major factor in its observed higher uptake than N2.

Possible relative motions between the gas molecules and the ions can be inferred from the density plots. Rotation of CO2 with the oxygen atoms moving toward the nearest anions, that is, inplane rotation following the arrow direction in Figure 7a, is predicted to be hindered because of the repulsion from the fluorine atoms bonded to P1. However, rotation of CO2 with the oxygen atoms into/out of this plane and toward the cations, as depicted in Figure 7c, seems to be favorable. These relative motions lead to localized P1 densities, whereas the P2 densities are more dispersed around the gas molecule. Furthermore, the P2 densities are slightly more localized closer to the oxygen atoms because of their preferred interactions. The average number of P1 and P2 atoms around a CO2 molecule within the 8 Å cutoff was calculated to be 2.6 for the P1 atoms in the anion and 3.1 for the P2 atoms in the cation (Figure 7d). The value of 2.6 is very close to 2.8 as predicted by the MD simulations for CO2−PF6 in the IL [bmim][PF6].18 The distances between the dense regions (shown as DAA′ and DBB′) can be compared to the RDF plots for P1−P1 and P2−P2 pairs in the OIPC. The separation DAA′ in Figure 7a is 8.2 Å, which is very close to the position of the first RDF peak of P1− P1 in Figure 8b. However, the distance DBB′ corresponds to the second peak in the RDF of P2−P2. Thus, the plots demonstrate the preferred coordination site of CO2 in the OIPC, consisting of nearest neighbors for anions, but second nearest neighbors for the cations. The distances are also indicative of the requisite separations to accommodate a CO2 molecule in the OIPC lattice. As discussed earlier, the structure (with vacuum interface) at 325 K is more disordered than the corresponding structure in the bulk solid (without interface). Unlike in the case of the disordered structure, the RDF peaks for P1−P1 and P2−P2 pairs in the ordered structure do not coincide with the separations DAA′ and DBB′. This can be seen in Figure 8a, wherein there is a wider spread of the P1−P1 distances about

4. CONCLUSIONS Classical MD simulations were used to study the uptake of CO2, N2, O2, and CH4 in the solid phases of the OIPC [P122i4][PF6]. The observed trend in gas uptake was CO2 > CH4 > O2 > N2. The extent of CO2 uptake was concluded to be determined by H


The Journal of Physical Chemistry B temperature and the relative amount of ordering in the solid structure, which in turn is expected to determine the free volume in the OIPC structure. The OIPC structure was, in general, relatively more disordered near vacuum and gas interfaces than in the bulk regions. The absorption and adsorption of CO2 were evaluated at four different temperatures. At 275 K, the solid exhibited the highest CO2 adsorption at the interfaces but negligible absorption within the solid because of a highly ordered OIPC structure. The highest bulk solubility was observed at 310 K, the OIPC phase (phase III) in which increasing dynamics and rotational disorder begin, and where the initial CO2 absorption creates even greater disorder after a long equilibration time. However, it appears that, when the disorder in the OIPC solid structure was already present before introducing the gas (i.e., at 325 and 365 K), the additional gas absorption that occurs does not further increase the disorder. This suggests that the availability of more free volume that results from a disordered OIPC structure favors CO2 absorption; however, as the material becomes more liquid like, the solubility in fact decreases because of an increase in temperature (which was also found for CO2 solubility in ILs). The CO2−OIPC structural analysis suggests that the average distances between the CO2 and the cation or anion is similar across all temperatures. These distances were concluded to determine the minimum prerequisites of intermolecular separation to accommodate CO2 in the OIPC lattice. When there is sufficient plasticity in the solid, the absorbed CO2 can create the volume required. One may then ask why the absorption was not greater at 325 K where there is even greater dynamics and presumably more free volume. The answer is likely that at higher temperatures, the adsorption of the CO2 (i.e., the surface excess) is less than that at 310 K; therefore, although there is sufficient dynamics and space available for CO2 to permeate into the OIPC at 325 K, the thermodynamic driving force is less. The simulations presented here suggest that the OIPC can be used as a potential membrane material for separating CO2 from N2, O2, and CH4 and that the behavior is phase-dependent. This has important implications in designing OIPC materials for optimum gas absorption.

■

■

ACKNOWLEDGMENTS

■

REFERENCES

The authors are grateful to the Australian Research Council for funding through the Australian Laureate Program FL110100013. This research was undertaken with the assistance of resources and services from the National Computational Infrastructure, which is supported by the Australian Government.

(1) Solomon, S. D. Q.; Manning, M.; Chen, Z.; Marquis, M.; Averyt, K. B.; Tignor, M.; Miller, H. L. IPCC: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change: Cambridge University Press: Cambridge, United Kingdom and New York, NY, USA, 2007. (2) Fuss, S.; Canadell, J. G.; Peters, G. P.; Tavoni, M.; Andrew, R. M.; Ciais, P.; Jackson, R. B.; Jones, C. D.; Kraxner, F.; Nakicenovic, N.; et al. Betting on Negative Emissions. Nat. Clim. Change 2014, 4, 850−853. (3) Kenarsari, S. D.; Yang, D.; Jiang, G.; Zhang, S.; Wang, J.; Russell, A. G.; Wei, Q.; Fan, M. Review of Recent Advances in Carbon Dioxide Separation and Capture. RSC Adv. 2013, 3, 22739−22773. (4) Baker, R. W.; Lokhandwala, K. Natural Gas Processing with Membranes: An Overview. Ind. Eng. Chem. Res. 2008, 47, 2109−2121. (5) Aaron, D.; Tsouris, C. Separation of CO2 from Flue Gas: A Review. Sep. Sci. Technol. 2005, 40, 321−348. (6) Ramdin, M.; de Loos, T. W.; Vlugt, T. J. H. State-of-the-Art of CO2 Capture with Ionic Liquids. Ind. Eng. Chem. Res. 2012, 51, 8149− 8177. (7) D’Alessandro, D. M.; Smit, B.; Long, J. R. Carbon Dioxide Capture: Prospects for New Materials. Angew. Chem., Int. Ed. 2010, 49, 6058−6082. (8) Yampolskii, Y. Polymeric Gas Separation Membranes. Macromolecules 2012, 45, 3298−3311. (9) Baker, R. W.; Low, B. T. Gas Separation Membrane Materials: A Perspective. Macromolecules 2014, 47, 6999−7013. (10) Bernardo, P.; Drioli, E.; Golemme, G. Membrane Gas Separation: A Review/State of the Art. Ind. Eng. Chem. Res. 2009, 48, 4638−4663. (11) Lin, H.; Freeman, B. D. Gas Solubility, Diffusivity and Permeability in Poly(Ethylene Oxide). J. Membr. Sci. 2004, 239, 105−117. (12) Bara, J. E.; Camper, D. E.; Gin, D. L.; Noble, R. D. RoomTemperature Ionic Liquids and Composite Materials: Platform Technologies for CO2 Capture. Acc. Chem. Res. 2010, 43, 152−159. (13) Blanchard, L. A.; Hancu, D.; Beckman, E. J.; Brennecke, J. F. Green processing using ionic liquids and CO2. Nature 1999, 399, 28− 29. (14) 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. (15) Lei, Z.; Dai, C.; Chen, B. Gas Solubility in Ionic Liquids. Chem. Rev. 2014, 114, 1289−1326. (16) Rogers, R. D. CHEMISTRY: Ionic Liquids−Solvents of the Future? Science 2003, 302, 792−793. (17) Giernoth, R. Task-Specific Ionic Liquids. Angew. Chem., Int. Ed. 2010, 49, 2834−2839. (18) 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. (19) 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. (20) Gu, Z.; Blanchard, L. A.; Hancu, D.; Beckman, E. J.; Brennecke, J. F. Proceedings of 5th International Symposium on Supercritical Fluids, Atlanta, USA, 2000.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.8b05444. Analysis of total gas uptake, CO2 density profile at 275 K, and spatial distribution of carbon atoms of the cations

■

Article

around N2 at 325 K (PDF)

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. Phone: +61392468383 (V.S.K.). *E-mail: [email protected]. Phone: +61392445114 (F.C.). ORCID

Vinay S. Kandagal: 0000-0003-0975-4725 Fangfang Chen: 0000-0002-8004-1720 Notes

The authors declare no competing financial interest. I


Article

The Journal of Physical Chemistry B (21) Kazarian, S. G.; Briscoe, B. J.; Welton, T. Combining ionic liquids and supercritical fluids: in situ ATR-IR study of CO2 dissolved in two ionic liquids at high pressures. Chem. Commun. 2000, 2047−2048. (22) Camper, D.; Bara, J.; Koval, C.; Noble, R. Bulk-Fluid Solubility and Membrane Feasibility of Rmim-Based Room-Temperature Ionic Liquids. Ind. Eng. Chem. Res. 2006, 45, 6279−6283. (23) Shannon, M. S.; Tedstone, J. M.; Danielsen, S. P. O.; Hindman, M. S.; Irvin, A. C.; Bara, J. E. Free Volume as the Basis of Gas Solubility and Selectivity in Imidazolium-Based Ionic Liquids. Ind. Eng. Chem. Res. 2012, 51, 5565−5576. (24) Carvalho, P. J.; Coutinho, J. A. P. On the Nonideality of CO2 Solutions in Ionic Liquids and Other Low Volatile Solvents. J. Phys. Chem. Lett. 2010, 1, 774−780. (25) Klähn, M.; Seduraman, A. What Determines CO2 Solubility in Ionic Liquids? A Molecular Simulation Study. J. Phys. Chem. B 2015, 119, 10066−10078. (26) Scovazzo, P. Determination of the Upper Limits, Benchmarks, and Critical Properties for Gas Separations Using Stabilized Room Temperature Ionic Liquid Membranes (SILMs) for the Purpose of Guiding Future Research. J. Membr. Sci. 2009, 343, 199−211. (27) Nguyen, P. T.; Wiesenauer, E. F.; Gin, D. L.; Noble, R. D. Effect of composition and nanostructure on CO2/N2 transport properties of supported alkyl-imidazolium block copolymer membranes. J. Membr. Sci. 2013, 430, 312−320. (28) McDonald, J. L.; MacFarlane, D. R.; Forsyth, M.; Pringle, J. M. A Novel Class of Gas Separation Membrane Based on Organic Ionic Plastic Crystals. Chem. Commun. 2016, 52, 12940−12943. (29) Jin, L.; Nairn, K. M.; Forsyth, C. M.; Seeber, A. J.; MacFarlane, D. R.; Howlett, P. C.; Forsyth, M.; Pringle, J. M. Structure and Transport Properties of a Plastic Crystal Ion Conductor: Diethyl(Methyl)(Isobutyl)Phosphonium Hexafluorophosphate. J. Am. Chem. Soc. 2012, 134, 9688−9697. (30) Chen, F.; Jin, L.; de Leeuw, S. W.; Pringle, J. M.; Forsyth, M. Atomistic Simulation of Structure and Dynamics of the Plastic Crystal Diethyl(Methyl)(Isobutyl)Phosphonium Hexafluorophosphate. J. Chem. Phys. 2013, 138, 244503. (31) Kandagal, V. S.; Chen, F.; Jónsson, E.; Pringle, J. M.; Forsyth, M. Molecular simulation study of CO2 and N2 absorption in a phosphonium based organic ionic plastic crystal. J. Chem. Phys. 2017, 147, 124703. (32) Iranipour, N.; Gunzelmann, D. J.; Seeber, A.; Vongsvivut, J.; Doherty, C.; Ponzio, F.; O’Dell, L. A.; Hollenkamp, A. F.; Forsyth, M.; Howlett, P. C. Ionic Transport through a Composite Structure of NEthyl-N-Methylpyrrolidinium Tetrafluoroborate Organic Ionic Plastic Crystals Reinforced with Polymer Nanofibres. J. Mater. Chem. A 2015, 3, 6038−6052. (33) Kamberaj, H.; Low, R. J.; Neal, M. P. Time Reversible and Symplectic Integrators for Molecular Dynamics Simulations of Rigid Molecules. J. Chem. Phys. 2005, 122, 224114. (34) Potoff, J. J.; Siepmann, J. I. Vapor−liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen. AIChE J. 2001, 47, 1676−1682. (35) Perez-Blanco, M. E.; Maginn, E. J. Molecular Dynamics Simulations of CO2at an Ionic Liquid Interface: Adsorption, Ordering, and Interfacial Crossing. J. Phys. Chem. B 2010, 114, 11827−11837. (36) Budhathoki, S.; Shah, J. K.; Maginn, E. J. Molecular Simulation Study of the Performance of Supported Ionic Liquid Phase Materials for the Separation of Carbon Dioxide from Methane and Hydrogen. Ind. Eng. Chem. Res. 2017, 56, 6775−6784. (37) Vujić, B.; Lyubartsev, A. P. Transferable Force-Field for Modelling of CO2, N2, O2 and Ar in All Silica and Na+ Exchanged Zeolites. Modell. Simul. Mater. Sci. Eng. 2016, 24, 045002. (38) Sun, Y.; Spellmeyer, D.; Pearlman, D. A.; Kollman, P. Simulation of the Solvation Free Energies for Methane, Ethane, and Propane and Corresponding Amino Acid Dipeptides: A Critical Test of the BondPMF Correction, a New Set of Hydrocarbon Parameters, and the Gas Phase-Water Hydrophobicity Scale. J. Am. Chem. Soc. 1992, 114, 6798− 6801.

(39) Shi, W.; Maginn, E. J. Atomistic Simulation of the Absorption of Carbon Dioxide and Water in the Ionic Liquid 1-n-Hexyl-3methylimidazolium Bis(trifluoromethylsulfonyl)imide ([hmim][Tf2N]. J. Phys. Chem. B 2008, 112, 2045−2055. (40) Huang, X.; Margulis, C. J.; Li, Y.; Berne, B. J. Why Is the Partial Molar Volume of CO2So Small When Dissolved in a Room Temperature Ionic Liquid? Structure and Dynamics of CO2Dissolved in [Bmim+] [PF6-]. J. Am. Chem. Soc. 2005, 127, 17842−17851. (41) Hong, B.; Panagiotopoulos, A. Z. Atomistic simulation of CO2 solubility in poly(ethylene oxide) oligomers. Mol. Phys. 2014, 112, 1540−1547. (42) Liu, H.; Dai, S.; Jiang, D.-e. Solubility of Gases in a Common Ionic Liquid from Molecular Dynamics Based Free Energy Calculations. J. Phys. Chem. B 2014, 118, 2719−2725. (43) Chen, F.; de Leeuw, S. W.; Forsyth, M. Dynamic Heterogeneity and Ionic Conduction in an Organic Ionic Plastic Crystal and the Role of Vacancies. J. Phys. Chem. Lett. 2013, 4, 4085−4089. (44) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (45) Butt, H.-J.; Graf, K.; Kappl, M. Physics and Chemistry of Interfaces, 3rd ed.; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, Germany, 2013. (46) Dang, L. X.; Wick, C. D. Anion Effects on Interfacial Absorption of Gases in Ionic Liquids. A Molecular Dynamics Study. J. Phys. Chem. B 2011, 115, 6964−6970. (47) Yan, T.; Li, S.; Jiang, W.; Gao, X.; Xiang, B.; Voth, G. A. Structure of the Liquid−Vacuum Interface of Room-Temperature Ionic Liquids: A Molecular Dynamics Study. J. Phys. Chem. B 2006, 110, 1800−1806. (48) Finotello, A.; Bara, J. E.; Camper, D.; Noble, R. D. RoomTemperature Ionic Liquids: Temperature Dependence of Gas Solubility Selectivity. Ind. Eng. Chem. Res. 2008, 47, 3453−3459. (49) Sato, Y.; Fujiwara, K.; Takikawa, T.; Sumarno; Takishima, S.; Masuoka, H. Solubilities and Diffusion Coefficients of Carbon Dioxide and Nitrogen in Polypropylene, High-Density Polyethylene, and Polystyrene under High Pressures and Temperatures. Fluid Phase Equilib. 1999, 162, 261−276. (50) Sato, Y.; Takikawa, T.; Takishima, S.; Masuoka, H. Solubilities and Diffusion Coefficients of Carbon Dioxide in Poly(Vinyl Acetate) and Polystyrene. J. Supercrit. Fluids 2001, 19, 187−198. (51) 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]− CO2Mixture. J. Phys. Chem. B 2007, 111, 4477−4487. (52) Zhang, Z.; Duan, Z. An Optimized Molecular Potential for Carbon Dioxide. J. Chem. Phys. 2005, 122, 214507.

J


Atomistic Simulation of Gas Uptake and Interface-Induced Disordering

Recommend Documents