Adsorption on Edge-Functionalized Bilayer Graphene Nanoribbons

Oct 8, 2012 - A combination of ab initio and classical Monte Carlo simulations is used to investigate the effects of functional groups on methane bind...
0 downloads 0 Views 4MB Size
Article pubs.acs.org/JPCC

Adsorption on Edge-Functionalized Bilayer Graphene Nanoribbons: Assessing the Role of Functional Groups in Methane Uptake Vinay S. Kandagal, Amardeep Pathak, K. G. Ayappa,* and Sudeep N. Punnathanam* Department of Chemical Engineering, Indian Institute of Science, Bangalore-560012, India ABSTRACT: A combination of ab initio and classical Monte Carlo simulations is used to investigate the effects of functional groups on methane binding. Using Møller−Plesset (MP2) calculations, we obtain the binding energies for benzene functionalized with NH2, OH, CH3, COOH, and H2PO3 and identify the methane binding sites. In all cases, the preferred binding sites are located above the benzene plane in the vicinity of the benzene carbon atom attached to the functional group. Functional groups enhance methane binding relative to benzene (−6.39 kJ/mol), with the largest enhancement observed for H2PO3 (−8.37 kJ/mol) followed by COOH and CH3 (−7.77 kJ/mol). Adsorption isotherms are obtained for edge-functionalized bilayer graphene nanoribbons using grand canonical Monte Carlo simulations with a five-site methane model. Adsorbed excess and heats of adsorption for pressures up to 40 bar and 298 K are obtained with functional group concentrations ranging from 3.125 to 6.25 mol % for graphene edges functionalized with OH, NH2, and COOH. The functional groups are found to act as preferred adsorption sites, and in the case of COOH the local methane density in the vicinity of the functional group is found to exceed that of bare graphene. The largest enhancement of 44.5% in the methane excess adsorbed is observed for COOH-functionalized nanoribbons when compared to H terminated ribbons. The corresponding enhancements for OH- and NH2-functionalized ribbons are 10.5% and 3.7%, respectively. The excess adsorption across functional groups reflects the trends observed in the binding energies from MP2 calculations. Our study reveals that specific site functionalization can have a significant effect on the local adsorption characteristics and can be used as a design strategy to tailor materials with enhanced methane storage capacity.

1. INTRODUCTION Development of an alternative, renewable, carbon-free, and sustainable energy source is one of the major challenges in the 21st century. An alternate energy source should aim at alleviating global warming and consequently climate change.1 In the absence of a hydrogen-based economy, transportation relies heavily on fossil fuel-based sources. In this scenario, the challenge lies both in reducing CO2 emissions, as well as other sources of pollution that commonly occur in diesel and petrolbased transportation vehicles. In this regard, natural gas (NG) has emerged as a viable transportation fuel in many parts of the world. Because of the large methane content, NG is a cleaner burning fuel and hence environmentally more friendly than gasoline. There are abundant NG reserves, spread across the world, and a large market is envisaged with NG, to meet the world’s growing energy demands. However, replacement of gasoline with methane involves overcoming many technical and engineering challenges. Methane is supercritical at room temperature (Tc = 190.3 K; Pc = 45.96 bar), and as compared to the high density liquid state of gasoline, a relatively large volume of NG is required to generate the energy equivalent of gasoline.2 Hence, various storage methods are employed with the primary aim of increasing the volumetric energy density. These involve storage by liquefaction (at 112 K and atmospheric pressure3), compression (at 200 bar), and storage by adsorption in a suitable adsorbent. The resulting fluid is © 2012 American Chemical Society

referred to as liquified natural gas (LNG), compressed natural gas (CNG), and adsorbed natural gas (ANG), respectively. The drawback with LNG is the high cost involved due to the associated cryogenic technology. Although CNG has been commercialized in different parts of the world, there still exist issues related to multistage compression requirements and safety. Among these storage methods, adsorption has the potential to become the most cost-effective and safe method.3 For adsorptive natural gas storage to achieve the volumetric energy density equal to that of CNG, the U.S. Department of Energy has set a storage density target of 180 V/V at 35 bar and 298 K. Microporous materials with their pore sizes in the nanometer to subnanometer regime and large surface areas are ideal candidates for achieving methane storage targets. In this regard, activated carbons, zeolites, metal organic frameworks (MOFs), and more recently covalent organic frameworks (COFs) are being pursued as potential materials for ANG technologies. Activated carbons are attractive due to their relatively lower costs and preparation routes from a variety of raw materials. They can attain high surface areas (greater than 3000 m2 g−1), high packing fractions, and can be prepared from less expensive Received: July 16, 2012 Revised: October 6, 2012 Published: October 8, 2012 23394

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403

The Journal of Physical Chemistry C

Article

have been reported in the study by Wood et al.11 for comparison with binding energies for graphene nanoribbons obtained from DFT theory. Functionalized benzenes offer both a computationally tractable system and a representative fragment of binding sites for situations encountered in carbons as well as MOFs. The MP2 calculations serve as a guide to choosing appropriate functional groups for adsorption isotherm simulations. In the second part of this Article, we carry out GCMC simulations to obtain adsorption isotherms and binding energies of methane on bilayer graphene nanoribbons with edges functionalized with OH, NH2, and COOH. A detailed comparison of commonly used classical forcefields with available CCSD binding energies for methane−benzene and methane−phenol is carried out to choose an appropriate forcefield for the GCMC simulations. Classical simulations include the effect of temperature and pressure on methane adsorption and quantify the enhancement in adsorption as the functional group composition is varied.

precursors. Thus, they fulfill most of the requirements for a good adsorbent. Activated carbons have highly disordered structures made up of stacks of broken graphene sheets with functional groups predominantly on the edges,4 which affects the surface chemistry of the adsorbent. These functional groups are formed during the activation process, and the specific chemistry depends on the activation method.5 The influence of surface functional groups on activated carbons has been known to influence the adsorption properties.6 Greinke et al.7 patented a method of enhancing methane uptake by halogenating activated carbon and report about a 21% increase in the methane uptake. Functional groups are also found in adsorbents such as zeolites and metal organic frameworks and can considerably influence the adsorption properties.8,9 Organic framework materials allow for a systematic manipulation of lattice structures and offer the possibility to tailor functional groups in the organic linkers. This strategy of selective functionalization and creation of additional adsorption sites has recently been identified as a potential method for improving adsorption capacity.9 A recent study by Wilmer et al.,10 which involved large-scale screening of potential MOF structures, indicates that methyl-, ethyl-, and propyl-containing functional groups improved the methane adsorption capacity when compared to amines and other halogen-based functional groups. Although it has been known that functional groups can influence the adsorption capacity of methane, the specific effect of the functional groups toward enhancement of adsorption has not been studied previously. As compared to experiments, computational studies are suited to quantifying the effect of functional groups on adsorption of gases. Recently, Wood et al.11 and Dutta12 used density functional theory (DFT) to compute methane binding energies for a wide variety of functional groups. The computations were performed with appropriately modified functionals to capture the dispersion interactions. Their results indicate that COOH, H2PO3, and CH3 groups enhance the binding energy when compared to hydrogen-terminated ribbon edges. Such studies that determine the binding energies around functional groups are expected to yield microscopic information on the affinity of binding sites toward methane and thus aid in the rational design of potential adsorbent structures. In this study, we use a combination of ab initio and classical Grand Canonical Monte Carlo (GCMC) simulations to investigate the effects of functional groups on adsorption of methane. Because the adsorption of methane on adsorbents such as activated carbon, MOF’s, and COF’s occurs through physisorption, accurate ab initio methods that capture the weak dispersion interactions between atoms are required. In this regard, the Møller−Plesset (MP) perturbation methods and the coupled cluster single, double, and triple excitation (CCSD) methods are the most popular.13−17 The CCSD-based methods, although more accurate, are computationally demanding and have not been as widely used, except for small systems.15 The MPn, where the n refers to the level of perturbation, are computationally less intensive and hence have been widely used to study physisorbed systems. In this study, the MP2 perturbation method is used to compute binding energies of methane on functionalized benzenes. We study the effect of electron-donating groups such as NH2, OH, and CH3 as well as electron-withdrawing groups such as COOH, H2PO3 on the binding energies and identify the methane binding sites. Some preliminary MP2 calculations on functionalized benzenes

2. AB INITIO RESULTS 2.1. Computational Details. The binding energies are obtained for methane interacting with functionalized benzenes, phenol (OH), benzoic acid (COOH), aniline (NH2), toluene (CH3), benzophosphoric acid (H2PO3), and nitrobenzene (NO2) using ab initio calculations. The calculations were carried out using second-order perturbation theory (MP2) and aug-cc-pVDZ basis sets.18 These basis sets have been previously used in similar studies involving benzene−methane16,17 and phenol−methane complexes.17 We carried out geometry optimization of the methane−benzene and methane-functionalized benzene complexes using MP2 level theory with the augcc-pVDZ basis set. Basis set superposition error (BSSE correction) was treated using the counterpoise correction method.13,19 The contribution from the BSSE correction is typically in the range of 2−10 kJ/mol for the systems investigated in this study. Methane can approach the functionalized benzene molecule from different directions. Because a systematic study of all possible configurations is not practically feasible, we study the binding energies with the side (S) and top (T) configurations as illustrated in Figure 1. In the S approach, methane approaches the functionalized group from the side in the benzene plane. In the T approach, three distinct initial configurations are studied. The methane molecule is placed directly above the benzene ring center (TB) or directly above the benzene carbon (TC) attached to the functional group or above the atom of the functional group attached to

Figure 1. Different initial configurations investigated for the methanefunctionalized benzene interaction. Nomenclature: TB, above the center of the benzene ring; TC, position above the carbon attached to the functional group; TG, position above the atom of the functional group attached to the benzene ring; and S, side position in the plane of the functional group. 23395

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403

The Journal of Physical Chemistry C

Article

Table 1. Comparison of Binding Energies for Methane-Functionalized Benzenes Obtained from Different Initial Configurations, Using MP2/aug-cc-pVDZ R

ETB, kJ/mol

−H −OH −NH2 −CH3 −NO2 −H2PO3 −COOH

−6.39 (−6.3a) −6.53 (−6.7a)[3.79] −7.14[3.76] −7.47[3.77] −7.35[3.69] −8.37[3.62] 6.71[3.79]

ETC, kJ/mol −6.47[4.06] −7.55[3.70] −7.77[3.57] −6.33[4.00] −8.24[3.74] −7.77[3.48]

ETG, kJ/mol

ES, kJ/mol

−7.37[4.12] −6.98[4.82] −7.77[4.07] −6.33[4.58] −8.23[4.58] −6.59[3.96]

−3.571 −5.5[2.59] −3.21[3.60] −2.26[3.71] −4.49[4.77] −7.18 [3.96] −5.4 (−5.3a)[4.36]

a

Data taken from Ringer et al., 2006.17 The numbers given in square brackets are distances in angstroms. Distances are between the C of CH4 and the center of the benzene ring (TB) and C of benzene attached to the functional group (TC). For TG, the distance is between C of CH4 and the central atom of the functional group; O in OH, N in NO2 and NH2, C in COOH, C in CH3, and P in H2PO3. In the S configuration, the distance is between C of CH4 and the central atom of the functional group as defined above for the TG distance.

Figure 2. Optimized configurations for the methane-functionalized benzene systems from MP2/aug-cc-pVDZ calculations for the top initial configurations (Table 1). In each case, the side view (left) and top view (right) are illustrated. The strongest binding is observed for H2PO3. In all cases, the methane molecule is found to be shifted toward the carbon atom of the benzene ring attached to the functional group (R). Color Code: White - H, Light Orange - O, Pink - N, Green - P.

the benzene ring (TG). These four different initial configurations were used to find the global minimum. The Gaussian 09 software20 was used for all of the calculations reported here. In ab initio studies of Tsuzuki and co-workers,14 MPn and CCSD calculations were carried out for the methane−benzene system with different positions and orientations of methane placed above the benzene molecule. In the optimized configuration, the center of the methane molecule is located above the center of the benzene molecule with the a C−H bond directed toward the center of the benzene molecule. Hence, we have chosen this orientation as the initial configuration in the present calculations. Geometry optimization calculations were performed with this initial configuration using an MP2 level of theory with aug-cc-pVDZ basis sets. Isolated benzene and methane were separately optimized using the corresponding basis sets and level of theory. In the initial configuration, the carbon atom of the methane molecule is located at a distance of 4.8 Å from the benzene plane for all of the top configurations as illustrated in Figure 1. 2.2. Interaction of Methane with Functionalized Benzene. The binding energies for the methane-functionalized benzene systems for the different initial configurations (TB, TC, TG, and S) are given in Table 1. The final optimized configurations for systems that showed the strongest binding are illustrated in Figure 2. As a reference, we also obtained the binding energy for the methane−benzene system.

In all cases, the TB, TC, and TG configurations yield a higher binding energy than the S configuration, and all of the functional groups considered have the effect of increasing the binding energy with respect to the bare benzene molecule in the top approach. This can be attributed to the greater interaction with the π cloud for the top (TB, TC, and TG) configurations when compared to the S configuration. In the S configuration, only H2PO3 is found to enhance the binding energy relative to benzene. From the calculated energies (Table 1), the H2PO3 group shows the highest binding energy with methane in the TB configuration. From the binding energies in Table 1, there is no clear correlation between the initial configurations (TB, TC, or TG) and the lowest binding energy state observed. For COOH and NH2, the minima are observed for the TC configuration, and in the case of H2PO3 and NO2, the minima were observed for the TB configuration. In the case of CH3 and H2PO3, the differences between the various intitial configurations were within 0.3 kJ/mol. From the minimum energy configurations, we observe that methane is located in the vicinity of the carbon atom attached to the functional group with the methane tripod directed downward (at an angle) toward the benzene plane. In the case of the S configurations, the influence of the functional group becomes more important as the interaction with the benzene π cloud is reduced. From the optimized geometries, we observe that the methane is oriented with the tripod toward the functionalized group in the S configuration. Additionally, in the case of −OH and −COOH 23396

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403

The Journal of Physical Chemistry C

Article

interatom and many-body terms contributing to the methane− methane interaction. In the all-atom five-site Kollmann model, each atom of the molecule is treated as an interacting site. The site−site interactions include both Lennard-Jones and electrostatic interactions. The parameters for both of these models are given in Table 2.

groups, we observe that one of the C−H bonds on methane aligns along the −OH groups as shown in Figure 3. Similar

Table 2. Lennard-Jones Parameters and Charges for OneSite and Five-Site Methane Models one-site model21

atom

Figure 3. MP2-optimized configurations for methane approaching (a) phenol and (b) benzoic acid from the side (S). In both cases, the methane tripod is directed toward the functional group with the fourth methane CH bond directed along the OH bond of the functional group as illustrated with the dashed lines.

σ (Å) ε (K) q (e)

orientation of the C−H bonds for methane has been observed in density functional studies of the interactions of methane with functionalized graphene nanoribbons.11,12 To briefly summarize our results from the MP2 calculations, we observe that functional groups on benzene can enhance the binding energy by 1−2 kJ/mol when compared to benzene. The greatest enhancement of 1.98 kJ/mol is observed for H2PO3 followed by COOH and CH3 where the enhancement is 1.38 kJ/mol. Our results for the benzene and phenol in the TB configuration are in good agreement with the published MP2 results at the aug-cc-pVDZ level.17 However, in the case of phenol, we find that the TG initialization predicts a stronger binding site located as illustrated in Figure 2.

C H C H C H

3.73

five-site model22 3.4 2.65 55.055 7.901 −0.66 0.165

148 0.0

The interactions between methane and carbon atoms of the graphene layer are obtained from the work of Do and Do.23 The forcefield parameters for both the one-site and the five-site models of methane interacting with carbon atoms of graphene have been tuned to match the experimental data. To model the interactions of methane with the atoms of the functional groups, two forcefields were considered, the Universal Force Field (UFF)24 and the Optimized Parameters for Liquid Simulations (OPLS).25 Together with the two models of methane, this results in four sets of forcefields describing interactions between methane and the functional groups (Tables 3 and 4). To make a proper choice, the binding

3. GCMC SIMULATIONS 3.1. Forcefields. The ab initio calculations of binding energy presented in the previous section showed the influence of various functional groups toward adsorption of methane. However, ab initio calculations do not take into account the effect of temperature on adsorption. To supplement these calculations, we performed GCMC simulations of adsorption of methane on a bilayer stack of graphene nanoribbons. The edges of the graphene nanoribbons are terminated with either a hydrogen atom or an organic functional group. The GCMC simulations require a forcefield that describes interaction of methane molecule with other methane molecules, interaction of methane with nonedge carbon atoms of the graphene, and interaction of methane with various atoms of the functional groups. In the above description, the carbon atoms at the edge of the nanoribbons are considered to be a part of the functional group. Its interaction with the methane molecules is modeled similar to the interaction of the carbon atom of benzene with methane molecules. To accurately simulate the adsorption of methane, the forcefield parameters need to be chosen carefully such that their predictions accurately match experimental data or ab initio calculations. The methane molecule can be modeled using either the united atom one-site model or an all-atom five-site model. In this work, we compare the forcefields from the one-site model proposed by Siepmann and co-workers21 and the five-site model proposed by Kollman22 with the ab initio results. The one-site methane model treats the entire molecule as a single interacting site. The form of the interaction potential is given by the Lennard-Jones (LJ) potential whose parameters have been estimated to reproduce the experimental vapor−liquid phase equilibrium data. The one-site model is an effective pair potential between methane molecules, which implicitly includes

Table 3. Classical UFF Forcefield Parameters Used for Benzene and Phenol atom

σ (Å)

ε (K)

q (e)

C H O

3.4 2.5711 3.1181

29.13 22.12 34.7221

0.0 0.0 0.0

energies predicted by these classical forcefields were compared to the CCSD(T) binding energies for the methane−benzene and the methane−phenol complexes reported in the work by Ringer et al.17 To compute the binding energies as predicted by classical forcefields, the methane molecule is placed at different positions along the axis passing through the center of the benzene ring (similar to the configuration studied in the ab initio calculations) for both the benzene−methane and the phenol−methane interactions. At each position, the interaction energy of methane molecule with the benzene/phenol is computed. For the five-site model, the computations were performed for 1000 random orientations, and the minimum value was chosen. The results from these computations for the benzene−methane interactions are shown in Figure 4a and b for the phenol−methane interactions. From Figure 4a and b, we observe that the OPLS forcefield combined with the Kollman five-site methane model provides the best match with the CCSD(T) (complete basis set) values for both the benzene and the phenol systems. The minimum value of the potential energy obtained with the UFF one-site methane model is in good agreement with the CCSD(T) values for both the benzene and the phenol systems. However, binding energies deviate from the CCSD(T) on both sides of 23397

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403

The Journal of Physical Chemistry C

Article

Table 4. Classical OPLS Forcefield Parameters Used for Calculating the Interaction of Methane with Different Functionalized Benzenesa atom C6H6 C6H5OH

C6H5NH2

C6H5COOH

C_b H_b C_b H_b C_g O H C_b H_b C_g N H(NH2) C_b H_b C_g C(COOH) O(O) O(OH) H

σ (Å)

ε (K)

q (e)

3.75 0.0 3.75 0.0 3.75 3.07 0.0 3.75 0.0 3.55 3.3 0.0 3.75 0.0 3.75 3.75 2.96 3.00 0.0

55.3541 0.0 55.3541 0.0 52.838 85.547 0.0 55.3541 0.0 35.2253 85.5472 0.0 55.3541 0.0 52.838 52.838 105.838 85.547 0.0

0.0 0.0 0.0 0.0 0.265 −0.7 0.435 0.0 0.0 0.18 −0.9 0.36 0.0 0.0 0.0 0.52 −0.44 −0.53 0.45

Table 5. Comparison of Binding Energies for Different Functional Groups Investigated in Classical GCMC Simulationsa methane with

classical forcefield (kJ/mol)

MP2/aug-ccpVDZ (kJ/mol)

CCSD(T)/CBSb (kJ/mol)

C6H6 C6H5OH C6H5NH2 C6H5COOH

−5.9 −6.21 −6.12 −7.39

−6.39 −7.37 −7.55 −7.77

−5.8576 −6.15

a CBS refers to the complete basis set limit. The five-site Kollman model for methane along with the OPLS forcefield was used to obtain the classical forcefield data. bCCSD data from Ringer et al.17

functionalized benzene molecule while finding the minimum in the potential energy. At each point, the interaction energy of methane molecule with the functionalized benzene is performed over 1000 random orientations. The reported CCSD values for benzene and phenol are in excellent agreement with classical OPLS values, and this concurs with the trends observed around the minima for the data shown in Figure 4. The MP2 binding energies are in general about 1 kJ/ mol larger than the classical binding energies and CCSD values. It is known that MP2 values tend to overestimate the binding energies when compared to the CCSD results.26 This trend is observed in Table 5 for the interactions of methane with benzene and phenol complexes. One noticeable difference between the MP2 and classical values is the ordering between OH and NH2 systems. The classical forcefields indicate a higher binding for phenol, whereas the MP2 numbers indicate the opposite. Nevertheless, the differences are less than 0.2 kJ/mol, and more accurate CCSD calculations could reconcile these trends. Because of the close agreement between the CCSD(T) data with the classical five-site Kollman−OPLS forcefield predictions, all GCMC simulations reported in this Article were performed with this forcefield. 3.2. Molecular Model of Adsorbent. Because the main motive for the present work is to bring out the effect of different functional groups on adsorption, we adopt a simple model for the adsorbent (carbon), consisting of two stacked graphene nanoribbons in AB-stacking as shown in Figure 5. In such a stacking arrangement, a carbon atom in the top layer is at the center of a hexagon formed by atoms of the lower layer. The nanoribbons have a width (distance between the edge carbon atoms of a graphene layer) of 15.62 Å. The stacking distance between the two graphene sheets is 3.35 Å. The edges of the graphene nanoribbons are terminated with either a

a

The same parameters are used for the functional groups in the calculation of adsorption isotherms. C_b and H_b refer to the benzene carbon and hydrogen, respectively. C_g refers to the benzene carbon attached to the functional group.

the minimum. In all cases, we observe that a lower binding energy is obtained with the five-site model for methane when compared to the one-site model. We note that the OPLS forcefield used for benzene and phenol is united-atom, wherein the H atom does not have explicit representation. In Table 5, we compare the binding energies for the functionalized benzenes using the classical Kollman five-site model and the OPLS forcefield, with binding energies obtained from our MP2/aug-cc-pVDZ calculations and with published CCSD data. The CCSD values for benzene functionalized with groups other than −OH are unavailable in the literature. The classical binding energies were calculated by placing the functionalized benzene in a box with x,y,z dimensions of 11.9, 8, and 8 Å, respectively, with the benzene ring lying in the x−y plane. The minimum value of the potential energy of interaction (binding energy) was evaluated by placing the methane carbon in grid points separated by a distance of 0.05 Å. The methane molecule probes the region both above, below, and on the sides of the

Figure 4. Comparison of classical binding energies for (a) methane−benzene and (b) methane−phenol interactions against CCSD data of Ringer et al.17 Calculations with the OPLS and UFF forcefields combined with the single site and five-site methane models show that the five-site model with the OPLS parameters for benzene/phenol yields the best agreement with the CCSD results. 23398

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403

The Journal of Physical Chemistry C

Article

Figure 5. The stacked nanoribbon model used for the classical GCMC simulations, consisting of two graphene layers. The carbons atoms of the top and the bottom layers are colored blue and brown, respectively. The edges can be functionalized with different functional groups, with varying concentrations. Shown here is the case of −COOH functionalization.

Figure 6. Isotherm and heat of adsorption plots for the different functional groups at 298 K. In all cases, COOH is found to yield the greatest enhancement when compared to H-terminated nanoribbons. The differences between the different functional groups are greater as the concentration of the functional group is increased from 3.125 to 6.125 mol %. 23399

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403

The Journal of Physical Chemistry C

Article

methane surface excess. Thus, the methane-functionalized benzene interactions are most representative of the graphene nanoribbon system at low pressures and low functional group concentrations. This is consistent with the nearly similar values of the methane surface excess and heats of adsorption at 3.125 mol % concentration for both the NH2 and the OH functional groups. Although the methane surface excess for OH is marginally higher than NH2, the differences between the heats of adsorption lie within 0.2 kJ/mol. The trends in binding energies obtained from the OPLS forcefield are reflected in the isotherms, implying that it is important to obtain reliable forcefields to generate accurate isotherms. Figure 7 compares

hydrogen atom or an organic functional group. The number of organic functional groups varies according to the concentration at which the simulations are performed. The system is periodic in all three directions with the two graphene nanoribbons placed parallel to the x−y plane. The box height in the zdirection is chosen such that adsorption on the graphene surface is not affected by its periodic images above (or below). In this work, the dimensions of the simulation cell are Lx = 30.62 Å, Ly = 29.5 Å, and Lz = 100 Å. This corresponds to 192 carbon atoms per graphene nanoribbon. Adsorption isotherms are calculated at 298 K for varying concentrations of the functional groups. Computations were carried out at three different mole fractions of the functional groups: 3.125, 4.2, and 6.25 mol %. This gives 12, 16, and 24 functional groups in the simulation cell, respectively. Apart from carbon, the other elements typically found in activated carbons are oxygen, nitrogen, and hydrogen. These elements can exist in/as various functional groups. The concentration of these functional groups can be tailored by a variety of chemical methods,4 and the concentration of oxygen can be increased to about 14 wt %.27 The oxygen and nitrogen concentrations considered in the present work range from 8.1 to 15.5 wt % and represent values that can be obtained in the laboratory for oxygen. Although nitrogen is generally observed at lower concentrations, we investigate the 15 wt % nitrogen in the NH2 functionalization for the purpose of comparison. The structure and coordinates for the functional groups were obtained from the MP2 optimization calculations of the isolated functionalized benzene molecules. In all calculations, the graphene nanoribbons and functional groups are rigid. 3.3. Adsorption Isotherms. Becauase the methane storage targets specified by the U.S. department of energy are at 35 bar and 298 K, GCMC simulations were performed at 298 K for pressures ranging from 10 to 40 bar to obtain the total methane loading in the simulation cell. The excess adsorption was then calculated from the total loading as follows: Γexcess =

Figure 7. Percentage increase in surface excess for different functional groups with respect to bare graphene at a pressure of 40 bar. Data are shown at the highest functional group concentration of 6.25 mol %.

the percentage increase in surface excess with respect to bare graphene (H functionalization) for 3.125 and 6.25 mol % at 40 bar for all of the diffferent functional groups studied. Comparing the relative increase in the surface excess due to functionalization shows that COOH functionalization increases adsorption by 44.5% at 40 bar pressure and highest concentration (6.25 mol %), with respect to bare graphene as shown in Figure 7. At both low and high functional group concentrations, the enhancement with the OH-functionalized ribbons is approximately twice that observed with the NH2 ribbons. In addition, we also see cooperative effects of the functional groups in enhancing adsorption. At low functional group concentration, the enhancements in the adsorption are proportional to the respective binding energies (Table 5). However, at higher functional group concentration, the enhancements are much higher than the corresponding increase in the binding energies. This is also reflected in the increasing values of the heats of adsorption with increase in functional group concentration (Figure 6). 3.4. Density Distribution. Density distributions yield information about the preferred adsorption sites. Density distribution of methane along the x direction (refer Figure 5) after averaging over the y and z directions is plotted in Figure 8 as a function of the functional group concentrations for both 10 and 40 bar. We observe that while the density on the graphene surface converges to a constant value on either side, there is strong local variation in the density along the graphene edges. At 10 bar and a functional group concentration of 3.125%, a small enhancement is observed due to the functional groups over bare graphene for the OH and NH2 groups. The improvement with OH over NH2 is clearly seen in the vicinity of the

∫V (ρ(r) − ρb ) dr acc

A

(1)

where ρ(r) is the density of methane in the vicinity of the interface, ρb is the bulk density, and A is the surface area available for adsorption. Vacc represents the accessible volume available for the adsorbent, which is the total simulation box volume subtracted by the volume occupied by the two graphene sheets. In this study, we used values of 1543.65 Å3 for the volume occupied by the graphene sheets and 921.58 Å2 for A. The isotherms for different mole fractions are shown in Figure 6 along with the heat of adsorption. In all cases, COOHfunctionalized nanoribbons yield the highest adsorption as compared to the other functional groups, resulting in an increase of methane adsorption by 20−44% when compared to the H-terminated nanoribbons, across the entire pressure range and compositions studied. The differences between the surface excess and the heats of adsorption between OH and NH2 are small, with a slightly higher excess adsorption observed for OH. This small increase is more apparent as the functional group concentration is increased. We point out that the binding energies obtained from the ab initio calculations are for the functionalized benzenes and only provide an indicative measure of the expected adsorption trends. The small increase in the binding energy for OH (−6.21 kJ/mol) over NH2 (−6.1 kJ/ mol) (Table 5) is consistent with the trends observed in the 23400

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403

The Journal of Physical Chemistry C

Article

Figure 8. Comparison of methane density distributions along the x direction, across different functional groups concentrations at 10 and 40 bar. Parts (a) and (c) are at 3.125 mol %, while parts (b) and (d) are at 6.125 mol %. COOH groups show the strongest affinity to CH4 when compared to OH and NH2. The local density enhancement around the COOH groups exceeds the density on bare graphene for 6.125 mol % at both 10 and 40 bar pressures. Because the system is symmetric, the data have been averaged across the center line.

graphene nanoribbon model, we study the effect of both functionalization chemistry and composition on the methane uptake characteristics. In the first part of the study, MP2 calculations were carried out on benzenes functionalized with OH, NH2, CH3, COOH, and H2PO3, and functionalization was found to enhance the binding energies from 1 to 2 kJ/mol when compared to benzene (−6.39 kJ/mol). In all cases, methane binding in a direction normal to the benzene plane was favored over binding to the functional group along the benzene plane. The greatest enhancement was observed for the case of H2PO3 where the binding energy was −8.37 kJ/mol. A similar binding energy for of −7.77 kJ/mol was observed for COOH and CH3 with values of −7.55 kJ/mol for NH2 and −7.37 kJ/mol for OH. We assessed the accuracy of the available classical forcefields for the interaction of methane with benzene and phenol by comparing the UFF and OPLS forcefields with both one-site and five-site methane models with reported CCSD results for similar systems. The OPLS forcefield with the five-site Kollman methane model was found to accurately capture the binding energy as well as the variation around the binding minima when compared to the CCSD results for the methane−benzene and methane−phenol interactions. In comparison, the MP2 binding energies for these systems were higher than the reported CCSD results consistent with expected trends and reported differences between these two methods. GCMC simulations were carried out using the OPLS forcefield with the five-site methane models on functionalized (OH, NH2, COOH) as well as H-terminated graphene nanoribbon model at a temperature of 298 K and pressures up to 40 bar. Adsorption isotherms and heats of adsorption indicate that COOH has the greatest enhancement toward methane uptake. At a mole fraction of 6.25%, an increase in the

functional groups. At low pressures (10 bar) and all functional group concentrations, the methane density in the vicinity of the functional groups (OH and NH2) is lower than the density observed on bare graphene (Figure 8a). A significant departure in this adsorption trend is observed for the case of COOH functionalization. At the lower pressure of 10 bar and a functional group concentration of 6.25 mol %, the density around the COOH group exceeds the density on bare graphene (Figure 8b). This enhancement is observed for both low and high functional group concentrations at a pressure of 40 bar (Figure 8c and d), and at 6.25 mol % the density around the COOH functional groups exceeds that of bare graphene by about 25%. At a concentration of 6.25 mol %, the increase in adsorption around the COOH functional group also elevates (to a small extent) the methane density at the bare graphene surface. Thus, the enhancement to the surface excess is dominated by the increased adsorption around the COOH groups. We point out that this enhancement in local density was not observed in other functional groups investigated in this study. As a consequence, COOH is found to increase the surface excess of methane by 44.5% when compared to 10.5% and 3.9% for OH and NH2, respectively, at a functional group concentration of 6.25 mol %. These adsorption sites are illustrated in the 2D density plots shown in Figure 9 where the elevated adsorption around the functional groups as well as on the bare graphene surface for the COOH-functionalized ribbons is readily observed.

4. CONCLUSION We have carried out ab initio and classical GCMC simulations to assess the effect of different functional groups on methane adsorption. By considering an edge-functionalized bilayer 23401

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403

The Journal of Physical Chemistry C

Article

Figure 9. 2D (x−y) density maps obtained by averaging in the direction normal to the graphene surface (z). Data are illustrated at 298 K and 40 bar at 6.25 mol % for the different functional groups investigated. The adsorption sites due to the functional groups are clearly observed in all of the functionalized ribbons. In the case of COOH, the local density at the functional groups exceeds the density on the graphene nanoribbons. All units in the color bar are in kmol/m3. Note that the simulation cell is repeated along the centerline of the periodic box illustrated in Figure 5, resulting in a periodic image with the functional groups facing each other as shown above.



methane surface excess of 44.5% was observed for the COOH groups when compared to H-terminated graphene nanoribbons. In comparison, the increases for OH and NH2 were considerably smaller at 10.5% and 3.7%, respectively. These trends were similar to the trends in binding energies observed for the classical forcefields, indicating that it is important to use accurate forcefields to differentiate between different systems and obtain adsorption isotherms, which are compared to experimental data. The enhancement in the case of COOH is reflected in both the higher heats of adsorption as well as the local enhancement in density around nanoribbon edges where the functional groups are present. Our study indicates that functional groups can have a significant effect on the local gas adsorption characteristics as illustrated with the edge-functionalized graphene nanoribbon model used in this study. We expect the adsorption trends illustrated in this work to be generic and extendable to other materials such as MOFs or COFs where recent attempts indicate that specific functionalization can be incorporated into the material framework to enhance the methane storage characteristics.9,10 Finally, we note that our results have direct bearing for improving methane uptake in high surface area graphene oxide-based carbons, which have recently been synthesized and tested for methane storage.28

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (K.G.A.); sudeep@ chemeng.iisc.ernet.in (S.N.P.). Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS



REFERENCES

K.G.A. would like to acknowledge Bharat Petroleum Corp. Ltd. for a research grant, and we would like to thank Shobhana Narasimhan and Brandon Wood for useful discussions while this work was in progress.

(1) Mandal, T.; Gregory, D. Annu. Rep. Prog. Chem. 2009, 105, 21− 54. (2) Lozano-Castelló, D.; Alcañiz Monge, J.; De La Casa-Lillo, M.; Cazorla-Amorós, D.; Linares-Solano, A. Fuel 2002, 81, 1777−1803. (3) Menon, V.; Komarneni, S. J. Porous Mater. 1998, 5, 43−58. (4) Marsh, H.; Rodriguez-Reinoso, F. Activated Carbon, 1st ed.; Elsevier Science: Oxford, 2006. (5) Menéndez-Daz, J.; Martn-Gullón, I. Interface Sci. Technol. 2006, 7, 1−47. (6) Chang, H. T.; Furuya, E. G.; Miura, Y.; Noll, K. E. Water Sci. Technol. 2000, 161−166. 23402

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403

The Journal of Physical Chemistry C

Article

(7) Greinke, R.; Bretz, R.; Mullhaupt, J. Method for Storing Methane Using a Halogenating Agent Treated Activated Carbon. U.S. Patent 5,372,619, 1994. (8) Malla, P.; Komarneni, S. J. Porous Mater. 1995, 1, 55−67. (9) Deng, H.; Doonan, C.; Furukawa, H.; Ferreira, R.; Towne, J.; Knobler, C.; Wang, B.; Yaghi, O. Science 2010, 327, 846−850. (10) Wilmer, C. E.; Leaf, M.; Lee, C. Y.; Farha, O. K.; Hauser, B. G.; Hupp, J. T.; Snurr, R. Q. Nat. Chem. 2012, 4, 83−89. (11) Wood, B. C.; Bhide, S. Y.; Dutta, D.; Kandagal, V. S.; Pathak, A. D.; Punnathanam, S. N.; Ayappa, K. G.; Narasimhan, S. J. Chem. Phys. 2012, 137, 054702. (12) Dutta, D. Methane Storage in Activated Carbon Nanostructures: A Combined Density Functional and Monte Carlo Study. M.Sc. thesis, Indian Institute of Science, 2010. (13) Cammi, R.; Bonaccorsi, R.; Tomasi, J. Theor. Chem. Acc. 1985, 68, 271−283. (14) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. J. Am. Chem. Soc. 2000, 122, 3746−3753. (15) Sherrill, C.; Takatani, T.; Hohenstein, E. J. Phys. Chem. A 2009, 113, 10146−10159. (16) Shibasaki, K.; Fujii, A.; Mikami, N.; Tsuzuki, S. J. Phys. Chem. A 2006, 110, 4397−4404. (17) Ringer, A. L.; Figgs, M.; Sinnokrot, M.; Sherrill, C. J. Phys. Chem. A 2006, 110, 10822−10828. (18) Dunning, T. J. Chem. Phys. 1989, 90, 1007−1023. (19) Boys, S.; Bernardi, F. Mol. Phys. 1970, 19, 553−566. (20) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; et al. Gaussian 09, revision A.1; Gaussian, Inc.: Wallingford, CT, 2009. (21) Martin, M.; Siepmann, J. J. Phys. Chem. B 1998, 102, 2569− 2577. (22) Sun, Y.; Spellmeyer, D.; Pearlman, D. A.; Kollman, P. J. J. Am. Chem. Soc. 1992, 114, 6798−6801. (23) Do, D. D.; Do, H. D. J. Phys. Chem. B 2005, 109, 19288−19295. (24) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024−10035. (25) Jorgensen, W.; Madura, J.; Swenson, C. J. Am. Chem. Soc. 1984, 106, 6638−6646. (26) Riley, K. E.; Piton̆aḱ , M.; Jurec̆ka, P.; Hobza, P. Chem. Rev. 2010, 110, 5023−5063. (27) Chingombe, P.; Saha, B.; Wakeman, R. Carbon 2005, 43, 3132− 3143. (28) Srinivas, G.; Burress, J.; Yildirim, T. Energy Environ. Sci. 2012, 5, 6453−6459.

23403

dx.doi.org/10.1021/jp307039m | J. Phys. Chem. C 2012, 116, 23394−23403