Noble Gas Separation using PG-ESX (X = 1, 2, 3 ... - ACS Publications

Article pubs.acs.org/JPCC

Noble Gas Separation using PG-ESX (X = 1, 2, 3) Nanoporous TwoDimensional Polymers Anna M. Brockway† and Joshua Schrier* Department of Chemistry, Haverford College, Haverford, Pennsylvania 19041, United States S Supporting Information *

ABSTRACT: Using planewave pseudopotential density functional theory and classical molecular dynamics simulations, we investigate the transport of noble gases through a family of two-dimensional hydrocarbon polymer membranes that generalize the “porous graphene” (PG) material synthesized by Bieri et al. by insertion of (E)-stilbene (ES) groups. We find that density functional theory overestimates the barrier height and empirical dispersion corrections underestimate the barrier height, compared to reference MP2/cc-pVTZ calculations on PG. The barrier height for noble gas transport is greatly reduced from PG to PG-ES1, but additional increases in the size of the pore in PG-ES2 and PG-ES3 lead to an attractive potential well instead of a repulsive barrier. Using the computed potential energy surfaces, we compute pressure- and temperature-driven tunneling probabilities of He isotopes, and refit an improved classical force-field. Using classical molecular dynamics simulations, we find that PG-ES1 has an He permeance of 6 × 106 GPU, which is 90 times greater than that of PG, and demonstrate high selectivity for He versus CH4, Ar, and CO2. These results indicate that PG-ES1 is a promising membrane material for separating He from natural gas, and separating He isotopes by tunneling differences.

■

for separations.6,7,12,14,17,18 These structures have the advantage of inherent, well-defined, pore geometries which can be densely and regularly spaced on the two-dimensional surface. A growing family of these two-dimensional polymers have been produced experimentally. The earliest of these is two-dimensional polyphenylene (referred to by the authors as “porous graphene” or PG, shown in Figure 1a) synthesized by Bieri et al. via Ag surface-promoted aryl−aryl coupling.27 Other existing porous structures include graphdiyne synthesized by Li et al. using Cusurface promoted coupling,28 two-dimensional dimethyl-methylene-bridged triphenylamine (2D-HT) synthesized on Ag surfaces by Bieri et al.,29 and two-dimensional 1,4-benzenediboronic acid (2D-BDBA) synthesized by Ourdjini et al. using Ag surface catalyzed coupling.30 Recently, two-dimensional polyester31 and borazine (BPL-2(H))32 polymers have been synthesized. This demonstrates the growing ability to synthesize two-dimensional nanoporous materials such as the ones we study here. Bieri et al. synthesized PG by coupling functionalized molecular building blocks to create a two-dimensional network, and they suggest that a similar bottom-up approach can use different building blocks to tune the pore size and spacing.27 Replacing the biphenyl-like bonding pattern in PG (Figure 1a) with a larger molecule will expand the pore, and the use of a chemical species with pi-bonding will maintain a twodimensional geometry. Here, we investigate the effects of

INTRODUCTION The growing demand for helium in a variety of industrial and scientific applications, combined with ongoing production shortfalls, has led to shortages of this element.1,2 Additionally, 3 He is crucial for neutron detection, but the stockpile of this isotope is depleted.3,4 Existing methods for performing both the chemical separation of He from natural gas and the separation of He isotopes require energetically costly cryogenic techniques. In principle, membrane-based gas separation methods have a much lower thermodynamic cost,5 and recent theoretical work has proposed the use of graphene-like twodimensional membranes for the chemical6,7 and isotopic6,8−10 separation of He gas. Here, we study a family of new twodimensional polymers, and identify one material, PG-ES1, which has improved properties for He separation. Theoretical studies investigating the use of two-dimensional membranes to perform chemical,6,7,11−19 ionic,20,21 and isotopic6,8−10,22 separations have taken two approaches. The first approach is to study transport through small pores introduced into graphene sheets.9−11,13,15,16,19−23 (A closely related strategy for gas separation using small pores introduced into carbon nanotubes is the topic of two recent papers.24,25) While these studies provide a general understanding of the effects of pore size and chemical functionalization, and recent experiments are making progress in the use of oxidative etching to produce small size-selective pores,26 it is unclear how to produce the precise atomic arrangements reported in these works, or a dense array of many such pores needed for practical applications. The second approachtaken hereis to use intrinsically nanoporous two-dimensional polymer structures © 2012 American Chemical Society

Received: October 15, 2012 Revised: December 8, 2012 Published: December 11, 2012 393

dx.doi.org/10.1021/jp3101865 | J. Phys. Chem. C 2013, 117, 393−402

The Journal of Physical Chemistry C

Article

computational methodologies taken in the literature to describe these systems, and to validate the assumptions used in these calculations. Second, we perform a systematic examination of the potential energy surfaces (PES) of the noble gases with PG, PG-ES1, PG-ES2, and PG-ES3 to elucidate trends in the interactions. Third, we use the computed potential energy surface for PG-ES1 to study the tunneling of He isotopes through these materials under both pressure- and temperaturedriven conditions. Fourth, we use the potential energy calculations to assess the validity of existing classical force fields, and show that these must be modified for certain noblegas systems. Finally, we perform molecular dynamics simulations to compute the permeance of He, Ne, and Ar through PG-ES1, and compare to previous results for CO2, N2, and CH4, to establish the capability of PG-ES1 as a membrane for the separation of He from natural gas.

■

COMPUTATIONAL METHODS Electronic Structure Calculations. Planewave pseudopotential density functional theory calculations using normconserving Trollier−Martins pseudopotentials37 and the Perdew−Burke−Ernzerhof (PBE)38,39 generalized gradient exchange-correlation functional were performed with Abinit 6.4.1.40,41 For the quantum mechanical calculations, PG-ESX unit cells were constructed using known molecular bond lengths and angles. The aromatic ring geometries were based on benzene, leading to C−C and C−H bond lengths of 1.421 Å and 1.07 Å, respectively, with a C−C−C bond angle of 120°. (E)-Stilbene bonds were used for the added groups in the PGESX structures, leading to central CC bond lengths of 1.356 Å; C−C bonds linking the central bond to the aromatic rings were set to 1.473 Å; C−H bonds on the alkene were set to 1.087 Å. CC−H bond angles were set to 120°, while angles for the C−CC bonds off of the central bond were set to 125°.33−35 Periodic calculations used a 2 × 2 supercell and were sampled only at the Gamma point in the Brillouin zone. For the comparisons between the finite (shown in Figure S1 in the Supporting Information) and periodic PG systems, we used a planewave energy cutoff of 150 Ryd. To determine a sufficient planewave energy cutoff for the remaining potential energy calculations, we performed calculations at increments of 10 Ryd from 50 Ryd to 180 Ryd with the He atom at the center of the PG system (shown in Figure S2 in the Supporting Information). The potential energy barriers calculated at 50 Ryd and 150 Ryd differ by only 10−4 eV. To validate the use of a low cutoff for the other noble gases, we performed calculations at 50 Ryd and 120 Ryd for the interaction of a larger-pore structure with all six noble gas atoms studied. Differences in the potential energy barriers were on the order of 10−5−10−6 eV for He, Ar, Kr, and Xe, and on the order of 10−4 eV for Rn. However, the difference was 0.003 eV for Ne; this larger value arises due to the construction of the Ne pseudopotential. Consequently, we used a planewave energy cutoff of 50 Ryd for computing the potential energy surface for He, Ar, Kr, Xe, and Rn, and a planewave energy cutoff of 120 Ryd for Ne in the work discussed below. Empirical van der Waals dispersion corrections were computed according to the D2,42 D3,43 and D3(BJ)44 methods described by Grimme and co-workers, using the dftd3 code.45 Because Grimme’s formalism is evaluated in real space, initially an expanded 6 × 6 supercell was constructed for the evaluation of dispersion corrections for the periodic PG system. Further investigation showed that distant atoms make negligible

Figure 1. Structures of the two-dimensional polymers considered in this work: (a) PG; (b) PG-ES1; (c) PG-ES2; (d) PG-ES3. van der Waals radius overlaps (gray circles indicate the proximal hydrogen atom van der Waals radii): (e) PG; (f) PG-ES1; (g) PG-ES2; (h) PGES3.

replacing the biphenyl-like bonding pattern in PG with (E)stilbene-like groups along one, two, or three of the hexagonal directions, as shown in Figure 1b,c,d, to yield structures denoted PG-ES1, PG-ES2, and PG-ES3, respectively. While (E)-stilbene at room temperature undergoes out-of-plane torsional vibration about the central CC bond,33−35 methyl substitutions on the phenyl rings reduce the magnitude of this torsional vibration.33 At low temperatures, this torsional vibration is reduced further to yield planar structures. The connectivity of the 2D polymeric structures in Figure 1b,c,d will limit this torsional disorder, providing a nearly flat twodimensional polymer. Qualitatively, the repulsive interaction experienced by an atom passing through the pore diminishes as the pore becomes larger. This can be rationalized by considering the overlap between the van der Waals radii of the hydrogen atoms lining the pore (shown as filled gray circles) and the corresponding radii of the noble gas atoms in the center of the pore (shown as outlines) in Figure 1e,f,g,h.36 PG has the most overlap with the van der Waals spheres of all the noble gas atoms, resulting in the largest potential energy barrier heights. The magnitude of this overlap is reduced for PG-ES1, yielding lower repulsive barriers than for PG. In PG-ES2, the van der Waals radii of the hydrogen atoms overlap with Ar, Kr, Xe, and Rn, but not with He or Ne. Finally, in the case of PG-ES3, the hydrogen atom spheres do not overlap with any of the noble gas atoms, and one expects to observe attractive potential energy wells for all of the species. This simple van der Waals sphere model does not consider the attractive dispersion forces, and, as we describe below, these corrections result in more complicated behavior for the larger PG-ES2 and PG-ES3 structures. This paper addresses several issues concerning these materials. First, we perform a set of benchmarking studies on PG to demonstrate the equivalence of several different 394

dx.doi.org/10.1021/jp3101865 | J. Phys. Chem. C 2013, 117, 393−402

The Journal of Physical Chemistry C

Article

contributions to the dispersion interaction, so a 2 × 2 supercell was used to calculate dispersion corrections for the periodic PG-ES1, PG-ES2, and PG-ES3 systems. Iterative Hirshfeld (Hirshfeld-I)46 atomic partial charges were computed using the self-consistent valence electron density in conjunction with all-electron atomic charge densities generated using the HF96 atomic Hartree−Fock code,47 as described by Glor et al.48 Our purpose here is only to describe qualitative properties of the charge distribution, but we note that Hirshfeld-I point charges more accurately reproduce molecular electrostatic potentials than alternative strategies such as Natural Population Analysis (NPA)49 and are more transferable between different molecular conformations than direct electrostatic potential fitting schemes.50,51 Tunneling Calculations. The transmission probability, t[E], of a noble gas atom as a function of kinetic energy, E = mv2/2, through the pore can be calculated from the DFT and dispersion-corrected potential energy surfaces using the onedimensional finite-difference method, as described by Cedillo.52 For simplicity, we consider only one-dimensional transmission, and so only the one-dimensional component (normal to the nanoporous material plane) for the kinetic energy. The thermally weighted transmission probability as a function of temperature, pw,m[T], is determined by weighting the transmission probability with the 1D Maxwell−Boltzmann kinetic energy distribution, pw, m [T ] =

∫0

∞

⎛ m ⎞1/2 −mv 2 /2k T B dv t[mv 2 /2]⎜ ⎟ e ⎝ 2πkBT ⎠

the complementary error function Erfc[x]. This allows us to simplify pw,m[T] to ⎛ m ⎞1/2 vh 2 pw, m [T ] = ⎜ t[mv 2 /2]e−mv /2kBT dv ⎟ vl ⎝ 2πkBT ⎠ 1/2 ⎤ ⎡⎛ mv 2 ⎞ 1 + Erfc⎢⎜ h ⎟ ⎥ ⎢⎝ 2kBT ⎠ ⎥ 2 ⎦ ⎣

∫

(2)

which leaves only the integral to be evaluated numerically. The ratio of the pw,m[T] functions for the two isotopes gives us an estimate of their selectivity by tunneling processes. We also considered the thermally driven tunneling process previously discussed in refs 8 and 10. r=

(PC,4 /PC,3) (PH,4 /PH,3)

⎛ p (TH) ⎞⎛ p (TC) ⎞ 1 1 w,4 ⎟⎜ w,3 ⎟e−(Δ3 −Δ4)( kBTC − kBTH ) =⎜⎜ ⎟ ⎜ ⎟ ⎝ pw,3 (TH) ⎠⎝ pw,4 (TC) ⎠

(3)

(4)

where Pα,i is the partial pressure of isotope i in the compartment having temperature Tα, the pw,i[Tα] are computed from eq 2 for each isotope, and Δi is the vibrational zero-point energy at the top of the barrier (i.e., of the activated complex). To determine the Δi, we performed single-point calculations in which we moved the He atom in the plane of the pore, calculating the potential energy surface at five and three additional positions in the positive “x” and “y” directions, respectively. The calculated values were assigned to symmetric positions in the negative “x” and “y” directions. The resulting potential energy surface was used to obtain the vibrational eigenstates. While this method neglects the coupling to vibrations of the pore atoms, which would be included in a complete harmonic normal-mode analysis, it more accurately accounts for anharmonic effects. Based on our previous calculations of this quantity in other nanoporous systems,8 coupling between the He atom and pore vibrations makes only a very small contribution to the normal modes, and can be ignored.10 We then used a finite-difference calculation to compute the bound states of this potential using the masses of the two He isotopes. Classical Molecular Dynamics. Following ref 17, the molecular simulations were performed using LAMMPS 2011.11.09,53,54 using the same structures, parametrization, and charges (when applicable) as in that study. To summarize: Interactions among PG-ES1 atoms were treated using the AIREBO force field.55 Making the nanoporous material rigid artificially increases the potential energy barrier for gas passage in many cases,19,25 so all PG-ES1 atoms are allowed to move freely in our simulations. Interactions between the gas molecules, and between gas molecules and PG-ES1 atoms, were treated using a Lennard-Jones (6−12) potential, truncated and shifted at 15 Å. The Lennard-Jones parameters for the PGES1 atoms were modeled using the benzene parameters of Fileti et al.,56 to describe the completely sp2-hybridized hydrocarbon framework of PG and PG-ES1. The LennardJones (6−12) parameters for the Ne and Ar atoms are taken from the parametrization of Vrabec et al.,57 and the parameters for the He atom are from the parametrization of Hirschfelder et al.58 As we will show below, it is necessary to refit some of the noble gas−hydrogen interactions for some systems to

(1)

where m is the mass of an individual atom or molecule, kB is the Boltzmann constant, and T is the absolute temperature. Following the discussion by Hauser and Schwertfeger,9 we set regions where t[E] < 10−7 to zero; this occurs for velocity vl and below. Similarly, for very high values of vh and above, t[E] = 1 (as shown in Figure 2); in this region, the transmission behaves classically and can be computed analytically in terms of

Figure 2. Model chemistry comparison: (a) He-PG potential energy surface; (b) helium tunneling probability. Solid lines indicate nonperiodic calculations, while dotted lines indicate periodic calculations. 395

dx.doi.org/10.1021/jp3101865 | J. Phys. Chem. C 2013, 117, 393−402

The Journal of Physical Chemistry C

■

RESULTS AND DISCUSSION Previous computational studies of gas transmission used either discrete molecules6,9,10 or periodic systems7,11,12 as models for nanoporous materials; the former is better suited for quantum chemistry methods, while the latter is better suited for solid state electronic structure calculations. Until now, no direct comparison between these two models has been made. Because PG has the smallest pores of the nanoporous membranes under investigation, it provides the limiting case for differences between the two models. The 2 × 2 supercell of PG is shown in Figure 1a, and the fragment (the same as used in ref 6) is shown in Figure S1 in the Supporting Information. The geometries are identical in both cases. The black lines in Figure 2a show the PBE results for the two models; the solid curve shows the nonperiodic finite fragment, and the dashed curve shows the periodic system. Qualitatively they are identical, and they differ by only 3.8 meV in barrier height. Consequently, the He atom tunneling probability, shown in Figure 2b, is also shifted by about 3.8 meV. From this we conclude that the finite fragment and periodic system are functionally equivalent. The -D2,42 -D3,43 and -D3(BJ)44 dispersion corrections of Grimme and co-workers introduce pairwise empirical longrange r−6 correction terms that could create a difference between the potential energy barriers computed for finite and periodic systems. To examine the source of the overall correction, we decomposed the contributions of the different atoms to the total dispersion correction for He atom passage, shown in Figure 3. The PBE-D2 correction was used as an

reproduce the potential energy barriers. These parameters are summarized in Table S1 in the Supporting Information. To obtain a tighter bound on the methane permeance, the simulations from ref 17 were extended to 150 ns using both a united-atoms model of methane by Vrabec et al.,57 and also an atomistic model by Stassen.59 The rigid geometry of the latter model does not allow for CH4 distortion as it passes through the pore. While quantum chemistry calculations indicate a large distortion when transiting through the very small pores of PG,6 only 1−2° deviations from tetrahedrality occur when transiting through larger pores comparable to those studied here,19 which justifies the rigid approximation used here and in other recent molecular dynamics simulations of comparable pore sizes.17,24,25 The membrane permeance and selectivity were obtained from the molecular dynamics simulations using the Langmuir adsorption based model recently introduced in ref 17. A series of equilibrium molecular dynamics simulations are performed under NVT conditions, and the total number of gas molecules that cross the barrier in both directions during the simulation of duration Δt is given by ⎛ αP ⎞ ⎟ rtotΔt = 2k 0[S0]⎜ ⎝ 1 + αP ⎠

(5)

where k0 is an effective barrier passage rate (accounting for both diffusion along the surface and passage through the pore), [S0] is the concentration of possible surface sites for adsorption, and (αP)/(1 + αP) is the fractional surface coverage, where α is a species-, surface-, and temperature-dependent adsorption equilibrium constant and P is the pressure in the system. For the present case, it is simpler to consider the Henry’s law limit of the Langmuir isotherm, where αP ≪ 1, which simplifies eq 5 to rtotΔt = 2k 0[S0]αP

(6)

having a single effective parameter, k0[S0]α, which can be obtained by fitting to the equilibrium MD crossing data and pressures (shown in Figure 8) obtained from several independent simulations. As in ref 17, crossing events were counted using the criterion of Du et al.,15 whereby a molecule is said to have crossed the barrier if it starts on one side of the barrier and then crosses to the other side and remains there for at least 1 ps. After the 1 ps time limit, a recrossing back to the original side is counted as another crossing event. A particle is defined as being outside the pore when it is ≥1 Å left or right of the polymer’s plane, but the results are insensitive to the particular cutoff used. In the Henry’s law regime described by eq 6, the permeance, pA, of a species A is constant with respect to pressure (and pressure differences) across the membrane, and is given by

pA =

k 0[S0]α aΔt

Article

Figure 3. Contributions of atom rings to PBE-D2 dispersion corrections in the nonperiodic system. The black dashed line indicates the total dispersion correction. All other lines represent the contribution of a single ring of atoms, as shown in the inset figure. Dotted lines correspond to rings of hydrogen atoms, and solid lines correspond to rings of carbon atoms.

explicit test case, but the results for the other methods are qualitatively the same. The largest dispersion correction terms originate from the hydrogen (dashed red) and carbon (orange) atoms nearest to the center of the pore. The second closest ring of carbon atoms (yellow) also plays a significant role in the dispersion interaction. However, atoms further from the pore contribute less than 10−4 eV, and thus the presence or absence of more distant atoms will have no impact on the dispersion correction. This indicates that the inclusion of the empirical dispersion corrections does not accentuate the difference between the periodic and finite-fragment potential energy surfaces. As seen in Figure 2a, the difference between the finitefragment (solid) and periodic (dashed) potential curves for

(7)

where a is the area of the membrane used for the simulation. The selectivity, sA,B, describes a membrane’s ability to discriminate between two species A and B. The ideal selectivity is simply the ratio of the permeances of two species, sA,B = pA/ pB. In the Henry’s law regime, the permeance is independent of the feed and permeate pressures. 396

dx.doi.org/10.1021/jp3101865 | J. Phys. Chem. C 2013, 117, 393−402

The Journal of Physical Chemistry C

Article

Potential Energy Surfaces. Having established the validity of the planewave pseudopotential DFT method, we calculated the potential energy surfaces for the four nanoporous membranes, shown in Figure 4. These potential energy surfaces were calculated for the 3 × 3 supercell of PG and for 2 × 2 supercells of the remaining structures. To validate the use of 2 × 2 supercells, barriers were also calculated for the 3 × 3 supercell of PG-ES1; these were identical to the 2 × 2 results to within 3 meV. The highest potential energy barriers occur for PG, shown in Figure 4a, and the relative ordering of the potential energy barriers increases with increasing noble gas size, as discussed above. MP2/cc-pVTZ data from ref 6 is also shown for He and Ne, and we see that in both cases this is in good agreement with the dispersion-corrected PBE. The high potential energy barriers for Ar, Kr, Xe, and Rn make PG essentially impermeable to the passage of these atoms. A considerable reduction in the potential energy barrier heights occurs for PGES1, shown in Figure 4b, consistent with the increase in the size of the pore shown in Figure 1f. He and Ne have low barriers of 0.072 and 0.123 eV, respectively, according to the PBE calculations. Because of the close proximity of the hydrogen atoms surrounding the pore, the dispersion corrections, which are on the order of 0.032 and 0.059 eV for these two species, have a proportionally large effect on the barrier heights. The dispersion corrections for Ar, Kr, and Xe lower the barrier heights by 0.179, 0.235, and 0.322 eV, respectively. However, these larger corrections are not enough to offset the high barriers of 0.725, 1.16, and 2.01 eV for these species. Further enlarging the pore to create PG-ES2 yields a qualitative change in the potential energy surface, shown in Figure 4c. The van der Waals radii of He and Ne do not overlap with the pore atoms of this structure, while those of Rn, Xe, Kr, and Ar do overlap (Figure 1c). Consistent with this, we see that the PBE potential energy barriers for He and Ne are attractive; this is further stabilized by the attractive dispersion correction terms. The dispersion correction for He and Ne is of the same order of magnitude as the well depth predicted by PBE. Although the simple van der Waals overlap model predicts a slightly repulsive interaction for Ar, the PBE results indicate that it is attractive even at the pore center. We attribute this to the polarizability of Ar, which is neglected in the simple hardspheres model formulation. The empirical dispersion terms act to further stabilize its potential energy well, lowering the well depth by 0.105−0.124 eV. The corrections are more than 2 orders of magnitude greater in energy than the (uncorrected) PBE result. The potential energy barriers for Kr, Xe, and Rn are more complicated. PBE predicts a repulsive potential energy barrier, consistent with the van der Waals sphere model. However, the polarizability of these gases leads to very large attractive dispersion interactions. These corrections more than offset the repulsive barriers: the more polarizable atoms end up with deeper wells, due to the dominating dispersion interaction. However, the dispersion corrections change the shape of the potential energy surfaces for these three noble gases in qualitatively different ways. In the case of Xe and Rn, the global minima are displaced from the center of the pore. This is consistent with the behavior of the dispersion interactions for He with PG shown in Figure 3. Near the pore center, the dispersion correction between the hydrogen atoms and the large noble gases is “turned off” in Grimme’s scheme; at distances further from the pore, this

each of the dispersion-correction models is nearly indistinguishable at this scale, and the magnitude is essentially the same as the 3.8 meV the difference between the (uncorrected) PBE curves. Grimme’s empirical dispersion correction models assume that partial charges on the membrane atoms do not change significantly as the molecular fragments approach each other.60 To validate the use of these corrections, Hirshfeld-I calculations were used to determine the partial charges in the periodic and nonperiodic porous graphene systems, and the results are shown in Figure S3 in the Supporting Information. For both finite and periodic systems, partial charges on the carbon atoms changed by approximately 0.015 e− when an He atom was placed in the center of the pore. The largest change in partial charges was observed on the inner-ring hydrogen atoms; this change is 0.027 e− for both the periodic and nonperiodic systems. An He distance of 4.0 Å from the pore is sufficient to restore the hydrogen partial charges to their original amounts. Changes in partial charges induced by the introduction of the He atom are small, justifying the use of Grimme’s dispersion correction formalism. Moreover, the fact that the change in partial charges is small and only affects the first two rings of atoms around the pore is consistent with only local interactions between the He atom and the pore. Having established the functional equivalence of the finitefragment and periodic model systems, we next compare the DFT results to the previously reported MP2/cc-pVTZ He-PG PES,6 which is taken to be authoritative. The MP2 calculation (yellow line in Figure 2a) is intermediate between the uncorrected-PBE (black) and the various empirical dispersion corrections (red, green, blue). Thus PBE overestimates the barrier height, and the various dispersion corrections underestimate the barrier height compared to the MP2 calculations. The potential wells range from approximately −0.043 eV for the PBE-D3 correction to approximately −0.0094 eV for the PBE calculations, and the minima are located between approximately 2.05 Å from the pore for the periodic PBE-D2 correction and 2.67 Å from the pore for the PBE calculations. The minima of the other models lie between these extrema. Tunneling probabilities were determined from the potential energy barriers using the finite-difference method (Figure 2b). The relative tunneling probabilities follow the same pattern as the potential energy barriers: lower barriers allow He atoms to pass through the pore with higher probability at lower incident kinetic energies. The dispersion correction shifts the tunneling curves by roughly 0.065 eV, as shown in Figure 2b. The periodic systems have marginally higher tunneling probabilities than the nonperiodic systems: 50% transmission occurs at roughly 0.0035−0.004 eV lower in the periodic systems, which again follows the slightly lower barrier heights predicted in those systems. In summary: The finite-fragment and periodic treatments of the potential energy surface agree to an order of magnitude less than the discrepancies between our different model chemistries. PBE overestimates the potential energy barrier height, compared to MP2. The Grimme empirical dispersion corrections tend to underestimate the potential energy barrier, compared to MP2. Of these, PBE-D2 and PBE-D3 give comparable results, and PBE-D3(BJ) is closest to the MP2 result. Generally, we can assume that (uncorrected) PBE will overestimate the true barrier height and the dispersion corrected results will underestimate it, allowing us to place upper and lower bounds on the transmission properties. 397

dx.doi.org/10.1021/jp3101865 | J. Phys. Chem. C 2013, 117, 393−402

The Journal of Physical Chemistry C

Article

For the PG-ES3 structure (Figure 4d), all of the noble gases have a simple attractive well whose depth increases going down the periodic table. The magnitudes of the well depths are all less than those of the corresponding atoms in PG-ES2, which can be understood by examining the van der Waals radius diagram discussed in the Introduction. In PG-ES2 (Figure 1g), the van der Waals spheres of the noble gas atoms are as close as they can be to the hydrogen van der Waals spheres without incurring a large repulsive interaction. Since the dispersion interaction falls off as r−6, the greater distance in PG-ES3 (Figure 1h) will mean a smaller attractive contribution from the hydrogen atoms. Tunneling Separation of He Isotopes with PG-ES1. Encouraged by the small barrier heightcomparable to that of the nitrogen-substituted pores recently studied by Hauser et al. for He isotope separation9,10we examined the tunneling of He atoms through PG-ES1. Figure 5a shows the transmission

Figure 5. (a) Thermally weighted transmission probability of 3He atoms through PG-ES1 as a function of temperature (eq 2); (b) 3 He/4He transmission ratio.

of 3He as a function of temperature for the different model chemistries. The thermally weighted transmission results from the four methods differ by about an order of magnitude at room temperature, and rapidly diverge as one approaches lower temperatures, e.g., the 77 K (liquid nitrogen) temperature in this plot. Figure S4a in the Supporting Information shows lower temperature results in the range of 10−30 K, where the PBED2 and PBE-D3(BJ) results give values of pw,3 > 10−5, but PBED3 gives pw,3 ≈ 10−11 and PBE gives pw,3 ≈ 10−27. This wide range of variation illustrates the unreliability of using results with empirical dispersion corrections for low-temperature tunneling predictions. Therefore, we primarily concerned ourselves with phenomena at 77 K and above, where this is not an issue. Figure 5b shows the ratio of the thermally weighted transmission for the two He isotopes. As expected, this

Figure 4. Potential energies of noble gas atoms approaching porous graphene systems. Solid points points/lines indicate the PBE results; dashed lines indicate the dispersion-corrected results. Insets show the same data close to the origin to emphasize details. (a) PG; (b) PGES1; (c) PG-ES2; (d) PG-ES3.

interaction is “turned on”. However, the long-range interaction with nearby carbon atoms leads to a large contribution at the center. 398

dx.doi.org/10.1021/jp3101865 | J. Phys. Chem. C 2013, 117, 393−402

The Journal of Physical Chemistry C

Article

to validate the molecular dynamics force-field parameters taken from the literature. As discussed in Computational Methods, we used the Lorentz−Berthelot combination rules to determine the noble gas−hydrocarbon interactions. Figure 7 shows a

asymptotically approaches 1 at high temperature, where classical (“over the barrier”) processes dominate. Values greater than 1 indicate selective transfer of 3He. At 300 K, the transmission ratio is between 0.6% (PBE, PBE-D3(BJ)) and 0.8% (PBE-D3). At 77 K, the transmission ratio is between 6.2% (PBE) and 7.4% (PBE-D3). This is notable because, even at this low temperature, the probability of He transmission through PG-ES1 is greater than the transmission probability through PG at elevated temperatures.6 Figure S4b in the Supporting Information shows this ratio in the range between 10 and 30 K. While the PBE-D2 and PBE-D3(BJ) models give ratios of 2.7 and 3.1, respectively, at 10 K, PBE and PBE-D3 give ratios of 138 and 119, respectively, at the same temperature. The true value is somewhere between these extremes. This leads us to offer a note of caution concerning the low-temperature results of Hauser et al.9,10 which used the B97D functional; subtle differences in the potential energy surface resulting from the choice of the empirical dispersion correction model can lead to dramatically different predictions. Figure 6 shows the temperature-driven isotope separation, recently discussed in refs 8 and 10. Only the PBE result is shown here; the dispersion corrected results (shown in Figure S5 in the Supporting Information) are qualitatively the same and differ primarily in the 3He transmission probabilities. For this temperature range, and for TC < TH, r < 1, indicating the dominance of zero-point energy differences. For TC = 273 K and TH = 373 K, r = 0.96 according to all of the models. The isotope separation can be magnified by working at lower temperatures, e.g., with TC = 77 K and TH = 300 K, r = 0.7 according to all of the models. However, at low temperatures (not shown in Figure 6) the models deviate; e.g., TC = 20 K and TH = 77 K, PBE and PBE-D3(BJ) predict that tunneling effects will dominate and give r = 1.4 and r = 1.3, respectively, whereas PBE-D2 and PBE-D3 predict that the zero-point energy will dominate and give r = 0.50 and r = 0.48, respectively. Classical Separation of He from Natural Gas with PGES1. The potential energy barrier data from Figure 4 allows us

Figure 7. Comparison on the DFT potential energies (symbols) for noble gases interacting with PG-ES1 to the literature Lennard-Jones interactions (solid lines): (a) He; (b) Ne; (c) Ar. In the latter, the black dashed line shows the potential energy obtained using a parametrization where the Ar−H interactions are fitted to match the PBE-D3(BJ) results for PG-ES1 (Table S1 in the Supporting Information).

comparison of the DFT data (symbols) to the potential energy computed in this way (solid lines), for the interactions between He, Ne, and Ar with PG-ES1. While the He and Ne parametrizations agree well with the DFT results, the Ar parametrizations yield a barrier which is much too high. To fix this, without modifying the adsorption properties of Ar on the carbon surface, we refitted the Ar−H parameters to reproduce the PBE-D3(BJ) results. The resulting potential energy (black dashed line in Figure 7c) better agrees with the DFT results; the parameters are given in Table S1 in the Supporting Information. A similar comparison was performed for the He− PG interaction (Figure S6 in the Supporting Information), which also required a refitting of the He−H parameters for PG. Fitting eq 6 to the crossing simulation data shown in Figure 8 yields the PG-ES1 permeances, which we report in units of mega-GPU (1 MGPU = 106 GPU = 3.3 × 10−3 mol cm−2 s−1 bar−1). The He, Ne, and Ar permeances are 6, 2, and 0.02 MGPU, respectively. For comparison, the He permeance of PG is only 0.07 MGPU (see Figure S7 in the Supporting

Figure 6. Temperature-driven isotope separation factors (eq 4) of He isotopes through PG-ES1, computed using the PBE tunneling potential and zero-point energies, indicated by black solid contours; blue dashed contours labeled on the right axis indicate the minimum 3 He thermally weighted transmission probability. 399

dx.doi.org/10.1021/jp3101865 | J. Phys. Chem. C 2013, 117, 393−402

The Journal of Physical Chemistry C

Article

membranes have He/CH4 selectivities