Rapid screening of nanopore candidates in nanoporous single-layer graphene for selective separations using molecular visualization and interatomic potentials

By virtue of being one-atom thick, nanoporous single-layer graphene (SLG) and other 2D materials can lead to the highest permeance and are thus considered the ultimate membrane.¹ In addition, by etching the right size of nanopores, it is possible to achieve molecular size-sieving, yielding outstanding selectivity.² Since separation processes account for more than 10% of the world’s energy consumption, SLG has tremendous potential to improve industrial energy efficiency.³ However, the development of 2D membranes faces challenges that need to be addressed to bridge the gap between theory and practical applications.^2,4–6

In particular, a fundamental understanding of nanoscale and sub-nanoscale mass transport and the ability to identify selective pores are important for membrane design to realize SLG’s promise.⁷ Targeting selective nanopores and estimating their selectivity are critical to optimize SLG membrane performance.⁸ However, studying mass transport across SLG at the atomic scale remains challenging on both experimental and theoretical fronts.⁹ Simulations are important for understanding the mechanisms of molecular transport through nanopores and sub-nanopores. Ab initio methods such as density functional theory (DFT) allow for calculations of transitional state barriers and are thus powerful tools to deduce transport mechanisms across nanopores and estimate transport rates.² However, their high computational requirements significantly reduce their range of applications. Molecular dynamics (MD) simulations have been intensively used to study mass transport across nanopores, where researchers have examined both gas-phase transport and liquid-phase ionic/fluid transport and investigated the effects of electrostatic interactions, nanopore shape, and/or adsorption.^2,10–12 For instance, electrostatic interactions between molecules and the nanopore rim, nanopore shape, and adsorption can severely affect the gas-phase transport behavior across SLG.^1,13–16 Ionic transport is also dependent on both electrostatic interaction and nanopore shape.^1,13–19 However, the simulation time for MD is still substantial. In fact, the computational cost of MD grows exponentially as the energy barrier of molecular permeation increases, which is associated with lower translocation probability, and hence, longer simulation times are required to fully capture the permeation behavior. For a typical simulation box containing 200 CO₂ gas molecules at 2 × 10⁶ Pa and 300 K, the simulation time is about 0.6 h for a computer with eight processors in parallel. However, for energy barriers up to 10 kJ mol⁻¹, the simulation time can reach 6 h, and for energy barriers reaching 20 kJ mol⁻¹, it can last 30 days.²⁰ As demonstrated above, the study of mass transport of a variety of molecules and ions across the numerous possible nanopores with various structures and functional groups is seriously hindered by computational requirements and the simulation time of conventional routes. Nevertheless, studying a wide range of nanopores is critical to identify which ones may be selective for specific pairs of molecules and also to understand how the heterogeneous pore size distribution (PSD) in an actual membrane determines its transport properties.

Herein, we report two simple, rapid-screening methods to quickly identify selective nanopores without relying on intensive simulations, analogous to the methods developed by Gounaris et al. to screen zeolites for selectivity.^21,22 Such methods, which consider the shape of the molecule and pore, have not been previously applied to molecular transport across nanopores in graphene or other atomically thin materials, and only effective pore and molecule sizes have been considered.^2,14 With the 3D representation of van der Waals radii for atoms, which accounts for the shape and size of nanopores, it was possible to quickly screen C-passivated nanopores for both gas- and liquid-phase separation. Using this tool, selective C-passivated nanopores for water desalinization/purification were identified, and shape-dependent mass transport across a given pore size distribution was examined. Furthermore, by taking electrostatic interactions into account, we also demonstrate a second method to screen selective non-C passivated nanopores for gas transport that accounts for interactions between partial charges on the atoms and their Lennard-Jones potentials. These screening methods with negligible computational requirements can be combined with MD simulations to enable faster assessment of transport across nanopores, facilitating our understanding and the development and adoption of selective nanoporous structures for real-world applications.

Visualization was used to rapidly screen selective nanopores involving two steps. First, optimized molecule and pore structures were created using 3D visualization software to represent the molecules that are to traverse across the nanoporous SLG. This was accomplished using Avogadro (version 1.2.0), a free, open-source advanced molecular editor and visualizer, which offers an easy-to-use chemical builder as well as different platforms for visualization and analysis with a powerful plugin architecture (http://avogadro.openmolecules.net).²³ To create the nanopores, carbon atoms can simply be removed from the SLG lattice. It is also possible to represent passivated nanopores with different kinds of atoms, such as H and N, by substituting the carbon atoms with other atoms [Fig. 1(a)]. For consistency, nanopores were named as follows: pore-Y-N_X, where Y is the number of missing carbon atoms, N is the atom used for passivation, and X is the number of those atoms.

A wide variety of force fields can be utilized to optimize molecule and nanopore geometries. This study used the universal force field (UFF) in Avogadro, which is a polyvalent force field that has parameters for every atom of the periodic table with an atomic number lower than 103 and has been used in a range of other simulation studies.^24–27 In addition, UFF is a non-reactive force field and therefore prohibits large geometric rearrangements such as covalent bond creation or breaking, which is essential for geometric optimization of both nanopores and molecules.²⁸ After shape optimization using UFF in Avogadro, a 3D representation of the van der Waals (vdW) radius for all atoms was displayed [Fig. 1(a)]. The vdW radius is directly derived from the force field parameters and is a good representation of the actual size of atoms relevant to transport across pores because it represents the distance of the closest approach and can thus serve as a good approximation of the real shape and vacancy area inside a specific nanopore and the actual size of the studied gas and liquid molecules.^23,29

Second, we overlaid the vdW-represented molecule and nanopore to assess whether the nanopore is permeable to a molecule. The molecule was placed at the center of the nanopore to examine the overlap between molecules/nanopores, and different molecule orientations were tested to check if at least one configuration fits within the nanopore. We expect the molecule to go through the nanopore (without a significant energy barrier) if there is no overlap between the molecule and the atoms at the pore rim based on the vdW representation for at least one orientation [Fig. 1(b)]. However, if the molecule does not fit in the nanopore, the molecule should not permeate readily through the pore [Fig. 1(c)]. This method assumes that the permeation of molecules through nanopores relies purely on the van der Waals size and provides a yes/no result as to whether a given molecule can readily permeate through a given nanopore.

To better characterize the size of molecules and nanopores, we introduced two shape metrics. The pore limiting diameter (PLD) was defined as the diameter of the biggest sphere that can fit inside the nanopore, taking into consideration the vdW radius [Fig. 2(a)]. In fact, PLD is a good representation of the permeable area of the nanopore for a molecule that just fits in the pore since it takes into account the jaggedness and irregularity of the nanopore rim. The molecule limiting diameter (MLD) was defined as the diameter of the circle that just encloses the smallest 2D projection of the gas molecule, also taking into account the vdW radius [Fig. 2(b)]. The ratio between PLD and MLD gives a good comparison between the permeable area of the nanopore and the size of the molecule.

We also report a second method that accounts for electrostatic and Lennard-Jones interactions. Unlike the 3D vdW visualization method, it accounts for both size and electrostatics. Energy barriers were quantified by calculating the contributions of van der Waals and electrostatics potentials. The potential barrier for a specific molecule–nanopore pair represents the energy barrier that the molecule needs to overcome in order to permeate through the nanopore, which is necessary for estimating the transmission probability across nanopores.² The overall potential encountered by a gas molecule when placed at the center of a specific nanopore is approximated by the sum of van der Waals potentials and electrostatic potentials of all atoms constituting the gas-phase molecules,¹⁵ 

where U_VDW (i, j) is the van der Waals potential between an atom of the gas-phase molecule (i) and a specific atom constituting the nanopore rim (j) and U_{electrostatic} (i, j) is the electrostatic potential between an atom of the gas-phase molecule (i) and a specific atom constituting the nanopore rim (j). The Lennard-Jones potential was used to approximate van der Waals interactions between each atom of the gas-phase molecule with atoms constituting the pore rim of nanopores,¹⁵ 

where ɛ₀ is the depth of the potential well, σ is the finite distance at which the inter-particle potential equals 0, and d_i,j is the distance between an atom of the gas-phase molecule (i) and a specific atom constituting the nanopore rim (j). The van der Waals parameters are listed in Tables S4 and S5 of the supplementary material.

To assess the electrostatic interactions between an atom constituting the gas molecule and an atom part of the nanopore rim, the gas-phase electrostatic interaction was assumed,³⁰ 

where q_i is the partial charge of one atom constituting the gas molecule, q_j is the partial charge of a specific atom part of the nanopore rim, and d_i,j is the inter-atomic distance between the two atoms. Partial charges of pore rims and gas-phase molecules provided by Sun et al. were used [Fig. 3(a)] and are listed in the supplementary material (Fig. S1 and Table S6).¹⁵ They were estimated by Sun et al. using density functional theory with the DMo13 module in Materials Studio.¹⁵ However, our approach is not restricted to performing density functional theory calculations for specific pores since any reasonable estimate of the partial charges and potentials may be used in Eqs. (1)–(3).

To calculate the center-to-center inter-atomic distances, the gas molecule was placed at the center of the nanopore, based on symmetry, using Avogadro [Fig. 3(b)]. In this study, only symmetrical nanopores and gas molecules were considered; for asymmetric pores or molecules, multiple configurations of the molecule in the pore may need to be calculated to identify the configuration that involves the lowest energy barrier. The distance between atoms in each atom-pair was calculated using ImageJ [Fig. 3(c)]. The orientation of the molecule was chosen such that it presents the lowest energy barrier to permeation. In the case of CO₂, the lowest energy barrier corresponds to the axis of the molecule being oriented perpendicular to the graphene, which was validated by calculating the energy barrier for a few different angles (for more complex geometries, it would be useful to qualitatively examine the repulsion and attraction between atoms to estimate the most probable orientation prior to interaction method analysis; a script can also be written to calculate the interaction energy across a range of orientations and determine the minimal energy pathway without having to perform full MD simulations). If the total potential was lower than the thermal energy at ambient temperature (k_BT), the gas molecule was assumed to permeate through the nanopore without any significant resistance. In this case, flux could be approximated by ideal gas flux,¹⁴ 

where ∆p is the driving pressure across the SLG membrane, M is the molar mass of the gas molecule, R_g is the universal gas constant, T is the temperature, A_membrane is the nanoporous SLG area, and A_p is the permeable area of the nanopores, defined as the vacant area in the 2D projection of the pore with the atoms being represented by their van der Waals sizes. Therefore, A_p/A_membrane represents the SLG membrane porosity. The permeance P is related to flux by J = PΔp.

However, if the total potential is higher than k_BT, the energy barrier would slow down the mass transport across the nanopore. To account for permeance or flux reduction, a correction factor taking into account energy barrier was used,² 

This transmission probability is based on the assumption that every molecule with a kinetic energy higher than the potential barrier can permeate through the nanopore. Molecular velocities are assumed to follow the Maxwell–Boltzmann velocity distribution.

The vdW-viz method offered the possibility of quickly identifying the smallest permeable (leaky) nanopore for a given molecule. This was achieved by progressively increasing the nanopore size until the molecule was able to permeate through the nanopore using vdW-viz. The nanopore size is described by its vacancy number, which is the number of missing C atoms in the pore. As the vacancy number increases, the number of plausible pore-isomers (i.e., pores with the same vacancy number but different structures) grows exponentially. For efficiency, only the two most probable isomers were considered for each vacancy number, as cataloged by Govind Rajan et al. based on their kinetic pathway analysis.³¹ If one of the pores is permeable to the studied molecule, the molecule was assumed to permeate through SLG at that vacancy number. If both nanopores are permeable, then the most probable pore was assumed to be the smallest leaky nanopore. The largest selective (non-leaky) pore that blocks the species to be retained (e.g., salt) while maximizing the permeation rate of smaller species (e.g., water) would generally be the pore with a vacancy number right below that of the smallest leaky pore. For practical membrane systems, the pore isomers at a vacancy number would likely exist as a distribution based on the kinetics of pore creation and thermodynamic equilibration (if any) such that if any pore isomer at that vacancy allows a molecule to permeate, the overall system at that vacancy would likely be leaky. In examples below, we assumed that pore occurrence probabilities follow the catalog of Govind Rajan et al.³¹ though the vdW-viz method can be applied to any distribution based on specific membrane fabrication methods.

In the case of identifying the smallest leaky nanopore for ion transport, the ions were assumed to be spherical, with their size described by their hydration shell diameter, while electrostatic interactions and dehydration were neglected.² Because of their large hydration shells, the smallest leaky nanopores could not be taken from the isomer catalog of Govind Rajan et al.; rather, they were guessed by progressively removing atoms from SLG until its pore-limiting diameter (PLD) was bigger than the cation hydration shell.³¹ 

To assess their accuracy, we employed the two screening methods to predict whether a specific nanopore is permeable to target gas molecules and used literature MD results as reference to compare predictions of the two methods, assuming MD studies to be representative of mass transport across SLG. Through this process, we could validate the methods and highlight the range of their applications.

The van der Waals 3D visualization method (vdW-viz) was compared against 13 different MD gas separation studies, in which the energy barrier and/or the permeance were reported, or where the permeance could be calculated using reported parameters.^{1,14–16,32–40} Of the 13 different MD studies, nine directly reported permeances or provided data from which permeance could be estimated, six studies reported energy barriers, of which five studies coupled DFT with MD to calculate energy barriers. A diverse set of nanopores with various shapes and passivation was studied to give a good representation of potential nanopores. 37 different nanopores with edges passivated by C, H, N, F, and/or O, and 11 different gas molecules were studied (H₂, N₂, CO₂, O₂, CH₄, H₂S, He, Ar, C₂H₄, C₂H₆, and CO), representing a total of 143 pore/molecule pairs (supplementary material, Tables S1–S3). Using the vdW-viz method, we predicted whether a target gas molecule can permeate through a given nanopore. The output of the screening method is binary, indicating whether a molecule can pass through a nanopore without a significant energy barrier.

Of the three gas molecule/C-passivated nanopore pairs where the energy barrier was reported, the vdW-viz method correctly predicted that the pore would be non-permeable to the molecule for the two cases where the energy barrier was higher than 0.1 eV and also correctly predicted that the pore would be permeable for the case where the energy barrier was lower than 0.1 eV (Table I and supplementary material, Table S2). At room temperature, the thermal energy k_BT is about 25 meV, and therefore, an energy barrier of 0.1 eV represents a reduction in the probability of transmission of a molecule across a pore by a factor of ∼0.02 and is therefore reasonable as a cutoff energy barrier. Most MD studies with C-passivated nanopores reported permeance or data that could be used to calculate permeance. Permeation predictions of the vdW-viz method were compared against permeances reported in the literature. The vdW-viz method successfully predicted that the pore would be permeable to the targeted molecule for the 55 cases where the normalized MD permeance by the ideal gas flux was higher than 0.01 (Table I). In addition, it also correctly predicted the non-permeability of the five pairs with normalized MD permeance to be lower than 0.01. Overall, of the 63 gas molecule/C-passivated nanopore pairs, we report 100% correct permeation predictions using vdW-viz (Table I and supplementary material, Table S2). These results show that the vdW-viz method is useful as a tool to predict the permeability of C-passivated nanopores to gas molecules.

To estimate the gas phase flux across nanopores, we plotted gas permeances of gas-molecule/C-passivated nanopore pairs reported by MD studies against the ratio of PLD/MLD, where each MD permeance is normalized by the ideal gas permeance [Fig. 4(a)]. We included only those MD studies that reported fluxes/permeances or where reported metrics allowed the proper determination of fluxes/permeances. Previous literature has mostly studied steric hindrance using the effective pore area/diameter, which, although taking into account vdW sizes, does not consider the nanopore/molecule shape that is relevant for selectivity.^2,14 The proposed metrics of PLD and MLD are thus more suitable for screening purposes. All permeable C-passivated nanopores reported in the MD studies had PLD/MLD >1, which is consistent with the vdW-viz method. Although gas permeance can be estimated using the ideal gas flux equation, the shapes of gas molecules and nanopores and gas molecules incident angles are important in molecular transport and, especially, in the steric regime where the energy barrier is not significant. As a first approximation to reduce the complexity of the problem, we approximated the steric coefficient, δ, which represents the probability of a molecule passing through a pore, given shape parameters, while neglecting adsorption and interactions as the remaining permeable area during the permeation of the gas pair molecule using a simple function of MLD and PLD,

The predictions matched the MD gas permeance data within an order of magnitude or better for PLD/MLD >1.5 [Fig. 4(a)]. For PLD/MLD <1.5, the permeance lies within an order of magnitude of the range defined by P_ideal and Eq. (6). Some of the discrepancies are explained by the simplification made in Eq. (6) that neglected the incident angle and the adsorption of gas molecules onto the SLG surface, which enhances the transport rate by offering another pathway of transport.¹⁴ As demonstrated by Sun et al., the surface flux can be higher than the direct flux, especially for gas molecules that can strongly adsorb onto the graphene surface such as N₂ and CH₄¹⁴ In addition, the surface flux becomes significantly higher for small PLD/MLD ratios because the surface flux is expected to scale as pore radius R_P and direct flux scales as R_P². It means that for smaller nanopores, the surface flux is more likely to dominate transport. Further work could be done on predicting gas adsorption impact on mass transport.

Assuming MD simulations to be representative of the true gas transport behavior, the result indicates that selective C-passivated nanopores for practical gas separation can be accurately predicted using the vdW-viz method. In addition, the results show that whether a gas molecule can pass through a C-passivated graphene pore is governed to a large extent by the vdW size, although electrostatics and surface adsorption can affect the overall transport rate.¹⁴ 

The vdW-viz method was also compared against 78 gas-molecule/non-C-passivated pore pairs reported in MD simulations (Table II and supplementary material, Table S3). Of the 34 gas molecule/non-C-passivated nanopore pairs where the energy barrier was higher than 0.1 eV, the vdW-viz method correctly predicted that the pore would be non-permeable to the molecule for 33 cases (Table II). Nevertheless, of the 22 molecule/non-C-passivated nanopore pairs where the normalized permeance is higher than 0.01, the vdW-viz method only predicted that the pore would be permeable for 11 cases. For example, in MD simulations, CO₂ is reported to go through the nanopore with a normalized permeance higher than 0.01; however, according to the vdW-viz method it should not go through.³⁹ In contrast, N₂ is not able to go through the same nanopore according to MD simulations;³⁹ however, vdW-viz predicts that it will go through. These inaccurate predictions suggest that parameters other than size are also essential for screening non-C-passivated nanopores. In fact, partial charges on rim atoms of non-C-terminated nanopores are larger in magnitude than in the case of C-passivated nanopores, and therefore, electrostatic interactions influence the permeation of molecules across non-C passivated nanopores [Fig. 4(b)].¹⁵ 

These results show that the vdW-viz method fails to be useful to predict whether a gas molecule will permeate across non-C-passivated nanopores. Successfully predicting whether a gas-molecule goes through a non-C passivated nanopore therefore requires considering electrostatics, which is discussed in Sec. III B.

To screen non-C-passivated nanopores, we utilized a second method that takes into account electrostatic interactions. The interaction method was checked against four different studies with four different nanopores and six different molecules, for a total of 20 pore/molecule pairs (supplementary material, Table S7).^15,38–40 Using the potential energy determined by Eq. (1) and kinetic gas flux Eq. (4), we could assess whether the targeted gas molecules can permeate through a given non-C-passivated nanopore. Only nanopores with N and/or H passivation were considered since other passivation could not be predicted from partial charges reported by Sun et al.¹⁵ The interaction method was compared against gas molecules/non-C-passivated nanopore pairs with reported permeance in the literature. Out of 20 cases, all except two permeation predictions matched with the MD simulations (supplementary material, Table S7). Sun and co-workers reported that H₂S could go through both pore-13–N₃–H₆ and pore-13–H₉ with permeances of 1.00 × 10⁻⁴ and 5.03 × 10⁻⁵ mol m⁻² s⁻¹ Pa⁻¹, respectively.¹⁵ However, according to our method, the energy barrier encountered by H₂S for those two nanopores are quite high (>0.1 eV), and the corresponding permeances were negligible (∼10⁻¹⁰ and 10⁻¹⁸ mol m⁻² s⁻¹ Pa⁻¹). The discrepancy can potentially be explained by inaccurate estimates of the potentials, or the very high H₂S adsorption onto the SLG surface in the MD simulation that could then impact H₂S permeation across the nanopores. Nevertheless, apart from H₂S, the interaction method successfully predicted gas permeances across non-C-passivated nanopores to within an order of magnitude (Fig. 5).

The gas interaction method described here is used to estimate a correction factor for ideal gas flux and is therefore not directly applicable to liquid-phase separation, which could include other effects beyond electrostatic interactions including dielectric exclusion, dehydration, and so on. For liquid-phase systems at the sub-nanometer length scale where few molecules can pass through the pore simultaneously, the solvent dielectric constant within the pore is not well-defined and would have to be adjusted based on more sophisticated simulations or experimental measurements before applying the electrostatic interaction method. Extending the interaction method to liquid-phase separation is left for future work.

In summary, we reported two methods to identify selective nanopores for molecular separation—vdW-viz method to screen selective C-passivated nanopores and interaction method to screen selective non-C-passivated nanopores. These two rapid screening methods are promising because of their good predictive ability as a screening tool and their simplicity.

By selective etching to control pore sizes and shapes, one can create tunable nanopores in SLG that can block micropollutants and salt while being highly permeable to water. Yet the vast range of possible nanopore shapes and sizes makes identifying selective nanopores particularly arduous. Given the good agreement of the predictions by the vdW-viz method with MD simulation results and its ability to explain the size-sieving phenomena, we employed the method to screen selective nanopores for water purification as a demonstration.

The vdW-viz method offers a way to quickly screen for nanopores that reject salt and micropollutants of interest while being permeable to water. This study examines four representative small organic pollutants: urea, a micropollutant coming from the use of pesticide in agriculture; N-nitrosodimethylamine (NDMA), a carcinogenic by-product formed during drinking water treatment; 1,4-dioxane, a carcinogenic by-product found at high concentration in wastewater from chemical factories; and ammonium perfluoro-2-propoxypropionate (GEN-X), used in the production of fluoropolymers, stain repellents, and as a component for aqueous film-forming foams [Fig. 6(b)].^41–44 The smallest leaky nanopore for each molecule was identified (Fig. 6). Knowing the smallest leaky nanopores for each species, it was possible to assess selective nanopores that are permeable to water while rejecting the selected micropollutants and salt (supplementary material, Fig. S2). These selective nanopores can be further studied through MD simulation and experiments.

This application demonstrates that using the vdW-viz method, selective nanopores could be screened for liquid-phase separations assuming size-sieving as the only separation mechanism. Another example of selective nanopore screening for aqueous rare-earth element separation is presented in the supplementary material, Sec. III B. Although these nanopores can be further studied through MD simulations and experiments, this method enables a faster and easier way to identify promising selective nanopores and therefore design 2D graphene membranes for targeted separations.

All fabrication processes for nanoporous graphene membranes will likely result in pore size distributions (PSDs) with pores of various shapes and sizes, which would impact mass transport. Assessing the effect of PSDs is challenging in MD simulations due to the wide variety of pores that are needed to be examined. In several studies that use simple models of transport across a pore size distribution, the shape of molecules and nanopores has been assumed as spherical and circular, respectively, and described by an equivalent diameter, thus neglecting the jaggedness and eccentricity of nanopores.⁷ However, slight changes in nanopore shape at the nanoscale can drastically affect macroscopic properties such as permeance.²⁰ Applying the vdW-viz method on experimental data enables us to understand how PSD affects transport while accounting for nanopore and molecule shape. We based our investigation on the work of O’Hern et al., where they studied the mass transport of potassium (K⁺) and Allura Red (AR), a larger non-spherical organic molecule on the order of 1 nm, across nanoporous SLG with different PSDs.⁷ We focused on the binary question on whether any of the pores present in the membrane under given pore creation conditions would allow species permeation. The nanopores were formed by gallium ion bombardment with subsequent oxidative etching by potassium permanganate, and different oxidative etching durations led to different PSDs.

To enable the application of the vdW-viz method, we first identified representative nanopores for different PSDs using the aberration-corrected scanning transmission electron microscopy (STEM) images of the SLG lattice obtained by O’Hern et al.⁷ Note that the STEM was performed on graphene transferred onto gold TEM grids, whereas the graphene for transport measurements was transferred onto polycarbonate track-etched membrane supports, and thus, the transfer and etching procedures might have resulted in slightly different nanopores.⁷ Nonetheless, the STEM images should still be representative of the actual PSDs present in the transport study. To obtain a rough estimation of the shape and size distribution of the PSD, a range of representative nanopores for each etching oxidative time reported by O’Hern et al. were represented using Avogadro (supplementary material, Sec. IV). Although the entire PSD could not be fully represented, a subset of the pores was imaged using STEM that spanned the pore size distribution in the PSD and was represented as an approximation of the PSD; 3–24 representative pores were used, depending on etch time. To estimate the nanopore shapes and molecular structure of representative nanopores, benzene rings were drawn around the STEM pore edges [Fig. 7(a)]. Subsequently, the nanopores were manually drawn in Avogadro [Fig. 7(b)]. The vdW-viz method was used to predict whether the molecules or ions could permeate through a given nanopore. Since hydrated K⁺ and (Cl⁻) ions could not be visualized in Avogadro, and given that the ions are approximately spherical, the permeability of a given pore to KCl was assessed by comparing the raw or hydrated ion diameters to the nanopore’s pore limiting diameter (PLD). Although our PSD representation is approximated by a subset of representative pores, by comparing the experimental results reported by O’Hern et al. with our estimated PSD, it was possible to expand further our understanding of mass transport across a given PSD. The calculation of the permeance for each PSD would require more sophisticated transport models and analysis of a greater number of pores in the pore size distribution, which can be a subject of future work.

AR transport across the PSD highlights the impact of nanopore shape on transport. Using the vdW-viz method, no AR leaky pores were identified (from among the representative pores) for oxidative etching times of 5 and 25 min. This matches perfectly with the experimental data since no increase in AR flux was reported for the same etching time.⁷ As a consequence, it means that AR is not able to permeate through nanopores smaller than its molecular size. In contrast, the mathematical model approximating the nanopores as equivalent circles with the same area that O’Hern et al. developed predicted AR to permeate through the PSD at 25 min.⁷ The discrepancy between the actual experimental/vdW-viz results vs the simplified model indicates that approximating nanopores as circles results in an overestimation of flux across PSD, especially for highly non-spherical molecules such as AR. Intuitively, it means that shapes of nanopores such as pore eccentricity and jaggedness would impact transport properties and thus affect separation. When emphasizing the nanopore and molecule shape using the vdW-viz method, it is possible to more accurately predict whether a given PSD is permeable to targeted molecules.

For an etching time of 60 min, some nanopores were predicted by vdW-viz to be permeable to AR. Although they represent the very tail of the pore size distribution, the AR flux across nanoporous SLG at 60 min was significant.⁷ This observation suggests that extreme nanopores dominate AR transport across nanoporous graphene. After 120 min of etching, the size of extreme pores did not change. Although the number of AR leaky pores increased, the AR flux reported by O’Hern et al. did not significantly change because mass transport was limited by the support resistance. The dominance of the tail of the PSD in governing AR transport can be explained by a simple scaling analysis that compares the transport resistance of graphene nanopores to that of pores in the underlying polycarbonate support (supplementary material, Sec. IV E).

In the case of KCl (where K⁺ and Cl⁻ ions have similar hydrated sizes), after 5 min of etching, none of the represented nanopores from the PSD were permeable to hydrated K⁺.⁷ Nevertheless, experimental results indicate an increase in KCl flux, which might be explained by the formation of K⁺ permeable nanopores that were not detected in the pore size distribution measurement. However, it is also possible that K⁺ (and Cl⁻) ions lose part of their hydration shell to permeate through smaller nanopores since the K⁺ raw cation diameter (1.33 Å) is much smaller than the hydrated K⁺ cation diameter (6.62 Å).⁴⁵ There would also be electrostatic interactions between the ions and nanopore rim; in the context of oxidative etching, negatively charged carboxylic groups are likely to be formed around the nanopores, as suggested by X-ray Photoelectron Spectroscopy (XPS).⁴⁶ Therefore, the favorable electrostatic interactions between the cation and carboxylate groups are likely to lower the energy barrier. The effective potential barrier can be expressed as the sum of the dehydration and electrostatic potentials, which is not considered in the current model.

For an oxidative etching time of 25 min, a few nanopores were found to be permeable to hydrated K⁺ cations.⁷ Similar to the case of AR, although they represented the tail of the total PSD, a significant KCl flux across SLG was reported. Similar to AR, extreme nanopores control ionic transport across PSD.

In summary, our analysis shows that the shape of nanopores is important to determine whether a PSD is permeable toward specific molecules and ions, and therefore, assuming nanopores as circular using the equivalent diameter/area would result in an overestimation of the permeability of the PSD since it does not take into account the nanopore geometry. By emphasizing the importance of size and shape using the vdW-viz method, it was possible to explain the experimental results reported by O’Hern et al. For small nanopores in the sub-nm scale, however, interactions such as electrostatic and dehydration may impact ionic transport, and these will need to be considered to determine the exact transport rates.

The theoretical study of mass transport across nanoporous graphene is limited by the high computational requirements of conventional simulation methods. In this work, we demonstrated alternative approaches to predict molecular selectivity through nanoporous graphene while accounting for nanopore and molecule shapes and interatomic interactions. First, using the 3D visualization of nanopores and molecules with their van der Waals radii, it was possible to assess whether a molecule can permeate through specific C-passivated nanopores, which was validated against gas-phase MD studies with 100% prediction accuracy over 63 molecule/pore pairs. Using this method, we could also quickly screen selective nanopores for water purification, desalination, and rare earth ion separation and also expand our understanding of mass transport across a given pore size distribution. The results suggest that extreme nanopores belonging to the tail of the distribution can dominate transport and highlight the importance of considering nanopore shapes instead of an “equivalent diameter” in determining the cut-off pore size. Accounting for electrostatic and van der Waals interactions, we also demonstrated another method to efficiently screen non-C-passivated nanopores. In addition, the gas-phase flux across non-C passivated nanopores could be successfully predicted to be within an order of magnitude. The reported methods fill the gap between simple estimates of permeation that treat pores and molecules as circles and spheres on the one hand and computation-intensive simulations on the other hand and could be extended beyond nanoporous single-layer graphene to other 2D membranes. Potential directions for future work include extending the interaction method to aqueous transport with electrostatic and hydration energies and considering the uncertainties in molecular and pore sizes due to effects such as conformational changes induced by molecule–pore interactions.⁴⁷ 

L.B. and C.-M.C. contributed equally to this work.

The supplementary material contains the list of MD simulations and their pore/molecule pairs, the parameters used in the interaction methods, the analysis of nanopores used to study the mass transport behavior across a given PSD, and the scaling analysis to estimate the contribution of the tail of PSD to transport.

L.B. acknowledges the Zeno Karl Schindler Foundation for the Masters Thesis Grant. C.-M.C. acknowledges support from the Croucher Foundation, the Abdul Latif Jameel Water and Food Systems Lab (J-WAFS) at the Massachusetts Institute of Technology, and the National Science Foundation Graduate Research Fellowship under Grant No. 1745302.

The data that support the findings of this study are available within the article and its supplementary material.