Efficient Kr/Xe separation from triangular g-C3N4 nanopores, a simulation study

Poly(triazine imide) or PTI is a promising material for separation of 85Kr/Xe mixture, thanks to its atom-thick nanoporous structure hosting triangular-shaped nanopores of ∼0.34 nm diameter.


Introduction
The development of materials to enable energy-efficient and straightforward processes can provide a solution to complex challenges in chemical separations. For example, an important separation challenge is purication of 85 Kr from a 10/90 mixture of 85 Kr/Xe, currently stored as hazardous waste. 1,2 Addressing this would reduce by ten times the storage volume of 85 Kr, and therefore reduce also associated risks. 3,4 Moreover, Xe recovered in a pure form or with negligible 85 Kr contamination can nd applications in ash and arc lamps, as propellants for ion thrusters in satellites, and as general anesthetics. [5][6][7][8] Currently, the separation of non-hazardous Kr/Xe mixtures is carried out predominantly by cryogenic distillation, relying on the relative volatility of liquid Kr and Xe. 7 When handling hazardous wastes from power plants, such as 85 Kr/Xe, reducing the process footprint is desirable, and membrane-based separation becomes very attractive. Membrane separation is one of the simplest industrial processes, requiring minimal maintenance. Moreover, by developing membranes yielding a high Kr permeance (dened as pressure-normalized ux), one can further reduce the membrane area needed. Nanoporous materials, hosting pores with a gap in the electron-density appropriate for sieving Kr (that has a kinetic diameter of 3.60 A) from Xe (kinetic diameter 3.96 A) can achieve high separation selectivity. [9][10][11][12] Recent progress in synthesizing atom-thick membranes, for example by using singlelayer graphene, has made it possible to realize large but selective gas permeance. [13][14][15][16] However, for effective sieving of Kr from Xe, two-dimensional nanoporous structures hosting uniform nanopores would be highly attractive. Such materials can be identied by screening databases of nanoporous materials. [17][18][19] A review of the literature indicates that carbon nitride materials can be very attractive to reach this goal. [20][21][22] Among them, g-C 3 N 4 is rapidly emerging for application in catalysis, 23,24 energy storage, 25 sensors, 26 and molecular separations. 27,28 The in-plane structure of g-C 3 N 4 hosts either a 1,3,5-triazine or a heptazine (1,3-,4,6,7,9,9bheptaazaphenalene) building unit. 29 In particular, poly(triazine imide) or PTI is quite attractive for gas separation because of its highly ordered nanoporous structure, composed of imide-bridged triazine units which form a triangular sub-nanometer-sized pore (Fig. 1). 27 A remarkable aspect of the mesh-like two-dimensional (2d) structure of PTI is its high pore density, 1.6 Â 10 14 pores per cm 2 , much larger than those typically realized in other 2d nanoporous membranes (such as nanoporous single-layer graphene), 30 making it promising for obtaining large permeances. The chemical and thermal robustness of PTI is also attractive for separation under harsh conditions. 31 Here, we investigate the performance of PTI nanopores in such useful membrane-based separation processes, employing vdW-DFT calculations, carefully benchmarked against exact exchange (EXX) plus RPA calculations, and molecular dynamics (MD) simulations. We indeed nd that PTI nanopores can be very attractive for Kr/Xe separation, with selectivity exceeding 10 000 at room temperature. This is accompanied by a high Kr permeance of 1000 gas permeation units or GPU (1 GPU ¼ 3.35 Â 10 À10 mol m À2 s À1 Pa À1 ).

DFT
To investigate the potential energy surface of Kr and Xe, vdW-DFT calculations are performed using the Quantum ESPRESSO distribution. 32,33 A 2 Â 2 supercell is used to ensure the decoupling of in-plane molecule-molecule interactions. Integration over the Brillouin zone is carried out with a uniform 6 Â 6 Â 1 k-point grid. A vacuum region of 40 A is employed in the z direction, to avoid interactions among periodic replicas. An energy cutoff of 60 Ry is used for the plane wave expansion of the wavefunctions. A kinetic energy cutoff of 480 Ry on the charge is used together with ultra-so pseudopotentials. 34,35 In order to carefully take into account the inuence of the van der Waals (vdW) interactions, ve different vdW approximations are tested: vdW-DF2, 36 Grimme-D2, 37 Grimme-D3, 38 the revised version of the Vydrov van Voorhis functional (rVV10), 39,40 and Tkatchenko-Scheffler (TS). 41 Since TS is not implemented in Quantum ESPRESSO when using ultraso pseudopotentials, calculations with this vdW approximation are carried out with norm-conserving pseudopotentials. 42 In this case, an energy cutoff of 100 Ry is used for the wave functions.

EXX/RPA
Using the implementation of the ACFDT formalism available in the Quantum ESPRESSO 33,43,44 distribution, we compute the total energies within exact-exchange plus the random-phase approximation. From a practical point of view, the main limitation of EXX/RPA is its computational cost, which is considerably higher than that of conventional LDA or GGA calculations. In order to use this computational-demanding method efficiently, some convergence tests are carried out (Fig. S2 †); the PES of gases is studied as a function of relaxation of the gas (Fig. S2a †), supercell size ( Fig. S2b †), k-point sampling (Fig. S2c †), pseudopotentials (Fig. S2d †), and cutoff values (Fig. S2e †). The PES is not affected signicantly when the supercell sizes and cutoff values are reduced, and when the pseudopotential is changed. Overall, it can be concluded that EXX/RPA calculations can be performed for the systems at hand at the G point using a PTI unit cell size and norm-conserving pseudopotentials.

Calculation of free energy prole
All classical MD simulations are carried out with the large-scale atomic/molecular massively parallel simulator (LAMMPS) package. 45 The interatomic interactions are described by the general AMBER force eld (GAFF). 46 Partial atomic charges on PTI are considered using the AM1-BCC method. 47 Kr and Xe are modelled as single center Lennard-Jones (LJ) particles. The LJ parameters used to consider non-bonded interactions, are tabulated in ESI Tables 1 and 2. † The particle-particle-particlemesh (PPPM) is implemented to consider long-range electrostatic interactions 48 and the vdW interactions are calculated within a cutoff distance of 12 A.
The dimensions of the simulation box in the x, y, and z directions are 25.9, 31.3, and 40.0 A, respectively. The PTI lattice is placed in the x-y plane and its atoms are allowed to move in all directions except the atoms on the external edges, which were kept xed to avoid any rigid shi of the PTI monolayer due to interactions with gases. In this way, it is ensured that the PTI sheet would not shi due to the interaction with Kr and Xe. Periodic boundary conditions are employed in all directions. For each simulation, the total energy of the system is minimized for 100 000 steps. Then, the system is equilibrated in a constant number, volume and temperature (NVT) ensemble for 1 ns with a timestep of 1 fs. A Nose-Hoover thermostat is employed to control the temperature of the system. 49,50 Then, in the restrained MD framework, umbrella sampling 51 is used to gather enough statistics in a reasonable period of time. To overcome the energy barrier and congurations which might be inaccessible in a reasonable amount of time, a harmonic restrain with the force constant of 5 kcal mol À1 A À1 is used every 0.5 A. Sampling is performed for 4 ns, and trajectories are collected every picosecond. Finally, the weighted histogram analysis method (WHAM) is used to unbias and combine the distribution of gases along the translocation trajectory function and construct the Helmholtz free energy prole. 52

Calculations of the concentration of occupied pores
In order to calculate the concentration of occupied pores, [C OX ], a simulation box with the dimensions of 25.9, 31.3, and 120.0 A Fig. 1 The structure of the nanoporous PTI lattice and its electron density isosurface (isovalue of 0.03 e A À3 ) revealing an electrondensity gap of $3.4 A. The dashed lines show the PTI unit cell.
in the x, y, and z directions, respectively, is dened (Fig. S3 †). We position the PTI lattice in the x-y plane at z ¼ 60 A, dividing the simulation box into two chambers (feed and permeate) of equal volume. In order to prevent the escape of the gas into the vacuum region due to the periodicity, a graphene layer is positioned at z ¼ 0 as a physical barrier. To accelerate the MD simulations and gather enough statistics in the 7 nssimulations, a gas pressure of 50 bar is used corresponding to 56 Kr or Xe atoms in the feed chamber while the permeate chamber is initially empty. In order to prevent vertical displacement of graphene and PTI, their coordinates are xed. The trajectories of gas atoms are collected every ten picoseconds to keep track of the number of adsorbed atoms. At the end of simulations, the average number of atoms in the adsorption region is obtained to calculate [C OX ] (Fig. S3 †).

Results and discussion
First, the PTI lattice is relaxed to its minimum energy conguration. The optimized lattice parameters are 8. 64, 8.65, 8.65, 8.64, and 8.67 A with vdW-DF2, Grimme-D2, TS, Grimme-D3, and rVV10, respectively, in good agreement with each other. To give a heuristic representation of the electron-density gap, the isosurface of PTI at a small value of 0.03 e A À3 is generated (Fig. 1).
An equilateral-triangle-shaped nanopore can be observed, with each side measuring 9.15 A. The nanopore can host a circle with a diameter of 3.40 A. This pore size is attractive for the separation of light gases based on their size, especially Kr from Xe by the size-sieving mechanism because the kinetic diameter of Kr is closer to the nanopore size. To investigate the adsorption behavior of gases on the porous PTI monolayer, a constrained relaxation is carried out. The gas is placed at several positions away from the sheet and the mass center of the gas is xed in the z direction (perpendicular to the sheet); the gas still has the freedom to relax in x and y directions to minimize the potential energy at a given distance from the PTI monolayer in the z direction. Since the potential energy surface (PES) for Kr and Xe are simple enough (Kr and Xe are just one single atom), it is not required to perform Nudge elastic band calculations. The convergence threshold on the total energy is 2 Â 10 À6 Ry and this criterion on forces is equal to 10 À4 Ry/a.u. The minimum PES for different gases passing through the PTI monolayer is explored by calculating the interactions energy with the PTI lm as: where E molecule+PTI is the total energy of the system when the gas interacts with PTI at different distances. E molecule and E PTI are the energy of the gas and the PTI monolayer, respectively, when they are isolated and do not interact. Fig. 2 shows the PES for Kr and Xe obtained using the ve considered vdW approximations. In order to make sure that the PTI monolayer does not shi rigidly due to the interactions with gases, the atoms on the edge of the PTI plane are kept xed. It is worth emphasizing that during these calculations all the PTI monolayer atoms (except the atoms on the edge of the plane) are free to move and relax. When the PTI lattice is treated as a rigid lattice, the PES prole did not change by a considerable amount (Fig. S1 †), indicating that the PTI lattice is reasonably rigid. The results indicate that the Kr adsorption energy in the middle of the pore (E ads , dened as the maximum attractive interaction energy when z s 0) is in the range of À0.06 to À0.17 eV, and the  32 Grimme-D2, 33 Grimme-D3, 34 rVV10, 35,36 and TS. 37 A height of 0 A corresponds to the plane of PTI.
adsorption height from the PTI plane is in the range of 1.60-2.00 A. Xe adsorption energy is in the range of À0.07 to À0.22 eV, and the adsorption height in the range of 1.80 to 2.50 A. The resulting activation barriers to translocate the nanopore (E act , dened as the difference in the interaction energy at z ¼ 0 and E ads ) for Kr and Xe are in the range of 0.14 to 0.22 eV, and 0.37 to 0.46 eV, respectively. In all cases, E act for Xe is higher than that for Kr ( Table 1). The discrepancy between the vdW approximations is more apparent in E ads than in E act .
To understand the effects of the vdW approximations on the sieving performance, the Arrhenius relationship for activation energies is used to calculate Kr/Xe selectivities, assuming that the translocation of the PTI nanopore would be the rate limiting step aer the gases are adsorbed on the PTI lattice. 53 The selectivity based on the Arrhenius equation is dened in eqn (2): where r is the rate of translocation and A Kr and A Xe are the prefactors for the Arrhenius relationship for Kr and Xe, respectively. E Kr and E Xe refer to the apparent activation energies (dened as E ads + E act ) for Kr and Xe, respectively; R is the universal gas constant and T is the temperature. Assuming that A Kr and A Xe are within the same order of magnitude, 54 an approximate gas pair selectivity within the temperature range of 250-400 K can be calculated (see Fig. 3).
Overall, all vdW approximations predict a high Kr/Xe separation selectivity, above 100 in all cases in the temperature range of 250-400 K (Fig. 3). Selectivities increase when the temperature is reduced, since Xe encounters a higher activation energy than Kr.
Despite the fact that all the vdW-DFT calculations give the same qualitative result, selectivities obtained with different functionals can differ quite dramatically, with roomtemperature values spanning a wide range that goes from 875 for Grimme-D2 to 7550 for vdW-DF2 (Fig. 3). Properly capturing vdW interactions with density-functional theory is complex and a relatively recent endeavor that still represents a challenging problem. 55 In order to have a higher-level description of van der Waals interactions, crucial to determine the gas-PTI monolayer interaction, and to have a qualitative estimation of the selectivity, for the rst time, we evaluated the energetics of the PTI complexes in the framework of adiabatic-connection uctuation dissipation theory (ACFDT), 56,57 allowing a better understanding of the performances of different vdW approximations for these systems. The EXX/RPA method has been proved to correctly capture vdW correlations, 58,59 and e.g., the results of physisorption binding energies between different types of molecules and surfaces show an excellent agreement with the values measured by experiments. [60][61][62][63][64][65][66] The RPA/EXX results reveal that the Grimme-D2, Grimme-D3 and TS approximations yield a good agreement PES (Fig. 4). The calculations show that vdW-DF2 overestimates the adsorption energy and on the other hand rVV10 underestimates the adsorption energy.  In absolute terms, the discrepancies with RPA appear to be larger for the two non-local functionals (rVV10 and vdW-DF2). This might suggest that the models for the effective polarizability of the two functionals, sharing an isotropic dependence on the charge density gradient, might give unprecise results in highly anisotropic congurations like gas-membrane interacting systems. However, further study is needed in order to elucidate this point. A minor role could also be played by the many-body interactions, fully included in RPA, partially included in Grimme-D3 through a three-body term and absent in the other approximations. 67 The transport of gases through two-dimensional nanopores can take place via direct gas phase transport or adsorbed-phase transport. 13,68 Attributing to the small nanopores of PTI, the gas transport is expected to proceed predominantly via the adsorbedphase transport mechanism, where gas rst adsorbs on the lattice, transport to the pore mouth in the adsorbed state, and nally arrive at the rate-limiting transition state at the center of the pore. Overall the gas pair selectivity is determined by the adsorption energy and relative size of the atoms/molecules with respect to the size of the nanopore. In the present case, based on the relative size of Kr (kinetic diameter of 3.6Å) and Xe (kinetic diameter of 3.96Å) with respect to PTI nanopore (3.4 nm), it is expected that size-sieving is the predominant separation mechanism. 69 However, while the vdW-DFT calculations reveal that PTI nanopores are promising for sieving Kr from Xe by the sizesieving mechanism, the selectivity calculations do not take entropic effects into account (loss of entropy of the gas in the adsorbed phase and in the transition state with respect to that in the gas phase). The calculations of entropic loss, for example using an enhanced method like umbrella sampling in the framework of classical MD simulations, allow one to accurately Fig. 5 The potential energy surfaces of Kr (a) and Xe (b) on a PTI monolayer, obtained with classical MD simulations to find out the best LJ parameters for Kr and Xe. The LJ parameters are trained in such a way that the potential energy surface shows a good agreement with the RPA results. The probability distribution of Kr (c) and Xe (d) along the translocation trajectories. Corresponding L (half-peak width of the probability distribution function) is shown. The Helmholtz free energy profile of Kr (e) and Xe (f) on a PTI lattice calculated with the umbrella sampling method at three different temperatures (300, 400 and 600 K) in the classical MD framework. predict the gas permeance using a gas transport model across 2d nanopores, as e.g., recently proposed by Strano and coworkers, 53,68 and more recently by Blankschtein and co-workers. 69 A unique advantage of this route is that one can avoid unaffordable long MD simulations to track the passage of gases. This is especially relevant when the activation barriers are high, making translocation of gases a rare event.
In this study, the LJ parameters are trained such that the resulting PES in the classical MD framework agrees with those from the RPA calculations ( Fig. 5a and b). Fig. 5e and f show the calculated Helmholtz free energy proles for Kr and Xe interacting with the PTI lattice at three different temperatures. The free energy barrier for both gases to translocate PTI nanopores increases with temperature. The barrier for Xe is larger than that for Kr, consistent with earlier observation from the vdW-DFT calculations.
In molecular transport across a nanoporous layer, it is wellestablished that when the interaction of a molecule at the center of the pore is repulsive, the corresponding transition state is the rate limiting step. Therefore, the translocation rate can be determined by measuring the rate of molecules crossing the transition state (k TST trans ). To calculate the rate of gases crossing the transition state, an Arrhenius-type equation can be used, and parameters of such a model can be obtained using the transition state theory: 70,71 where DE is the energy barrier which corresponds to the activation energy needed for the translocation of gases from the adsorbed state above the pore to the center of the pore or to the transition state (as illustrated in Fig. 5e). Based on Fig. 5a where L is the half-peak width of the probability distribution of gases along the translocation trajectory function ( Fig. 5c and d), M is the molecular weight of gases, and DS is the corresponding entropy changes for the pore translocation. DS can be calculated by taking the rst derivative of the Helmholtz free energy changes with respect to the temperature: where DF is the Helmholtz free energy changes when the gas moves from the adsorbed state to the transition state ( Fig. 5c  and d). Plotting DF as a function of temperature and by measuring the slope of the curve, the entropy loss is calculated (Fig. 6). The entropy losses are 38.0 and 39.7 J mol À1 K À1 for Kr and Xe, respectively. The entropy changes are referred to as loss because the degree of freedom of the activated complex at the transition state is lower. Interestingly, in the context of separation, these calculations point that the entropic losses are similar for Kr and Xe, perhaps attributing to the fact that both of them have a similar shape (spherical atoms). Since the entropic losses for Kr and Xe are comparable, their effect on the Kr/Xe selectivity is not dramatic ($15% at room temperature). Based on the entropic loss, we calculated the translocation prefactor, A trans , and the rate of gas translocation through the center of the pore, k TST trans , which are reported in Table 2 and discussed later. Once k TST trans is known, the gas ux can be obtained as: Fig. 6 The Helmholtz free energy changes of Kr (left) and Xe (right) at different temperatures when the gas arrives at the transition state from the adsorbed state.
where X, O, and O X represent the gas phase, unoccupied nanopore and occupied nanopores, respectively. The results show that on average, 4.5 Kr and 5.0 Xe atoms adsorb on top of the PTI lattice in the simulation box. Considering the pore density in PTI, the corresponding [C OX ] are equal to 9.6 Â 10 À7 and 1.1 Â 10 À6 mol m À2 for Kr and Xe, respectively (Table 2), based on which the ux and permeance are calculated. Finally, the Kr/Xe selectivity can be extracted by taking the ratio of the respective uxes.
The key parameters related to adsorption and translocation of Kr and Xe are summarized in Table 2. Briey, the concentration of adsorbed atoms corresponds to site occupancy of 37 and 41% for Kr and Xe, respectively. This sub-monolayer adsorption at a high pressure of 50 bar is expected for the noble gases. The adsorption of gas atoms sets the stage for the translocation event during which the adsorbed gas loses entropy. Since the pre-exponential factor, A TST trans , is an exponential function of entropy loss (eqn (3)), it decreases from the typical 10 13 s À1 to 2.0 Â 10 9 and 1.5 Â 10 9 s À1 for Kr and Xe, respectively. Despite this, thanks to the presence of a single transition state in the translocation of the gases across the atom-thick PTI nanosheets, a high Kr permeance of 1000 GPU at 300 K is obtained. Owing to the relatively high activation energy for Xe to translocate the PTI nanopore, Xe permeance is low (0.05 GPU), leading to a selectivity of 20 000.
Overall, the calculated Kr permeance (1000 GPU) and Kr/Xe selectivity (20 000) from PTI at room temperature is signicantly better than those from the state-of-the-art membranes demonstrated for Kr/Xe separation. Typically, selectivities in the range of 6-45 and Kr permeance in the range of 20-450 GPU have been reported. For example, 2-5 mm thick silicoaluminophosphate-34 (SAPO-34, chabazite zeolite) lms hosting a nominal pore size of 0.38 nm could separate Kr/Xe with a selectivity of 23 and Kr permeance up to 26 GPU. 1 Another study using 3-mm-thick SAPO-34 membranes yielded Kr permeance of 360 GPU and selectivity of 45 which is the best reported membrane performance to date. 9 The lower selectivity in the case of SAPO-34, with respect to that predicted from PTI membrane, can be attributed to relatively lower energy barrier of Xe to transport across 0.38 nm sized SAPO-34 pores compared to 0.34 nm sized PTI pores. Similarly, using another nanoporous material with a pore size of 0.38 nm (aluminophosphate-18 or AlPO-18), only a small Kr/Xe selectivity of 6.4 could be obtained. 10 The separation of Kr/Xe with nanoporous materials hosting smaller pores (0.34 nm) has also been attempted, involving zeolitic imidazolate framework-8 or ZIF-8 membranes. However, owing to the lattice exibility of ZIF-8, which tends to accommodate larger molecules in ZIF-8, a small selectivity of 14.2 was realized. 11 The extremely high Kr/Xe selectivity from PTI nanopore makes it especially attractive for the separation of 85 Kr from a 10/90 mixture of 85 Kr/Xe, targeting near pure 85 Kr and Xe streams using a PTI membrane. This will allow one to safely use the puried Xe for commercial applications. The recent development in obtaining single-layer PTI nanosheets 72 will allow one to fabricate a few nanometer thick lms, similar to what has been achieved in the case of graphene and metal-organic framework (MOF) nanosheet-based membranes. 73,74

Conclusions
This study reveals the high potential of PTI-based membrane, thanks to $0.34 nm-sized nanopores, for separation of Kr/Xe mixtures. In order to ensure accurate predictions on gas selectivity, which are critically inuenced by small changes in the adsorption and the activation energies, we performed an extensive comparison of several vdW-corrected functionals within DFT and benchmarked them against accurate, reference EXX/RPA calculations. We highlight that differences in commonly used vdW functionals can lead to dramatic changes in selectivity measurements, raising accuracy concerns for some of these that require further care of practitioners in the eld. In order to consider nite-temperature effects, the entropic contribution in the transition state has been determined from Helmholtz free energy calculations at different temperatures, obtained with restrained MD simulations where the force-eld parameters have been trained on the EXX/RPA potentialenergy prole. Finally, gas permeances are predicted by calculating the translocation coefficient from single-layer PTI using transition-state theory. We obtain, from our simulations, a Kr/ Xe separation factor exceeding 10 000 at room temperature, which would allow one to obtain nearly pure 85 Kr and Xe streams.
Overall, the method developed here, vdW-DFT calculations benchmarked against EXX/RPA calculations to predict the energy barrier and entropic penalties for gas translocation across a nanopore, can be applied to a diverse set of two-dimensional nanoporous lms, e.g., nanoporous graphene, C 2 N, polyheptazine imide, etc., for calculating the molecular separation performance. Apart from the high accuracy, the method used here to predict permeance is also attractive in comparison to the classical MD simulations when the energy barrier for the translocation event of at least one component is high enough to make the permeation a rare event, in which case, one has to conduct time-consuming and expensive long MD simulation trials.

Conflicts of interest
There are no conicts to declare.