Cutting Materials in Half: A Graph Theory Approach for Generating Crystal Surfaces and Its Prediction of 2D Zeolites.

Scientific interest in two-dimensional (2D) materials, ranging from graphene and other single layer materials to atomically thin crystals, is quickly increasing for a large variety of technological applications. While in silico design approaches have made a large impact in the study of 3D crystals, algorithms designed to discover atomically thin 2D materials from their parent 3D materials are by comparison more sparse. We hypothesize that determining how to cut a 3D material in half (i.e., which Miller surface is formed) by severing a minimal number of bonds or a minimal amount of total bond energy per unit area can yield insight into preferred crystal faces. We answer this question by implementing a graph theory technique to mathematically formalize the enumeration of minimum cut surfaces of crystals. While the algorithm is generally applicable to different classes of materials, we focus on zeolitic materials due to their diverse structural topology and because 2D zeolites have promising catalytic and separation performance compared to their 3D counterparts. We report here a simple descriptor based only on structural information that predicts whether a zeolite is likely to be synthesizable in the 2D form and correctly identifies the expressed surface in known layered 2D zeolites. The discovery of this descriptor allows us to highlight other zeolites that may also be synthesized in the 2D form that have not been experimentally realized yet. Finally, our method is general since the mathematical formalism can be applied to find the minimum cut surfaces of other crystallographic materials such as metal-organic frameworks, covalent-organic frameworks, zeolitic-imidazolate frameworks, metal oxides, etc.

Introduction

Two dimensional (2D) materials are quickly gaining attention as promising materials in a wide variety of applications.[1] While the pre-eminent 2D material (graphene) has a thickness of one atom,[2] more classes of 2D materials have been studied in recent times (i.e., van der Waals layered structures).[3] Atomically thin or layered crystals can also be considered 2D materials since they maintain long-range connectivity in the direction of two unit cell vectors but not the third. We hypothesize that answering the following question can yield insights into preferred crystal faces and hence the propensity of crystals to form 2D-like structures: how can a 3D crystal be cut into two separate partitions by severing a minimum number of bonds or, given a pairwise potential, by severing the minimum total bond energy per unit area? Such concepts have already been used to rationalize the crystal dissolution of a zeolite by atomic force microscopy,[4] and more recently it has been used in coarse grain crystal growth modeling.[5] Answering this question is, generally speaking, the analogous problem of minimal cuts in graph theory.[6] Given a connected graph, a minimum cut typically seeks to split the graph into two partitions by removing a minimum number of edges or, in the case of a weighted graph, removing edges whose total weight is minimal. Hence if a crystallographic material is interpreted as a graph with atoms for nodes and bonds for edges, we can use graph theory techniques to determine the surface termination for a given Miller plane that minimizes the number of cut bonds, i.e., the number of dangling bonds. As phrased here, this problem is not constrained to one particular class of materials. Initially, however, zeolites present the perfect class of crystals to apply this minimum cut formalism since the Si–O bonds have an extremely narrow range of distances and hence each bond can be thought of as approximately equivalent in strength. For the remainder of this work we focus our discussion on zeolite materials. However, we stress that the formalism of this graph theory approach could be applied to many other crystalline materials, including but not limited to metal oxides, zeolitic-imidazolate frameworks (ZIFs), and metal–organic frameworks (MOFs), by introducing edge weights that address the different bond strengths in more complex materials.

Zeolites are crystalline solids whose microporous void structures have made them ubiquitous in a variety of commercial applications, most notably as catalysts in the petrochemical industry and as adsorbents in small molecule separation processes. It is estimated that the global zeolite molecular sieve market will surpass USD 35 billion by 2024.[7] Just as the majority of experimental zeolite research has focused on the performance of bulk 3D materials in catalysis, separations, and other applications, so too has the majority of computational research focused on 3D materials, from hypothetical structure generation[8−10] (including successful prediction of novel materials[11]) and synthetic descriptors[12−17] to performance prediction and high-throughput screening.[18−21] However, the research on 2D zeolite materials has accelerated in recent times, especially over the past decade.[22−25] For example, 2D zeolites have demonstrated improved catalytic performance[26−28] and shown potentially greater separations efficiency[29−33] over their 3D counterparts. Both applications benefit from improved mass transfer in thin 2D zeolite materials. For catalysis, catalyst deactivation by coking, which is a major problem for 3D zeolite frameworks, can be significantly suppressed with 2D zeolite catalysts. For separations, the diffusion time for molecules to pass through the 2D zeolite molecular sieves can also be dramatically reduced. 2D zeolites have also been used as precursors to 3D zeolites for which no other synthetic procedure is known.[34] While the 2D forms of a relatively small number of International Zeolite Association (IZA)[35] zeolites have been discovered, their significant potential in a variety of applications merits further investigations to uncover novel 2D zeolite structures. In this work we show that in silico design of 2D zeolites can help further this goal.

Computational investigations of 2D zeolite structures have also intensified recently, mainly to understand the interactions between and the reassembly process of layered 2D precursors that are generated in the Assembly, Disassembly, Organization, Reassembly (ADOR) process.[36−39] We aim to build upon these computational studies to develop a high-throughput screening technique and a descriptor to understand whether an IZA zeolite is likely to demonstrate a stable 2D form. To do this we take advantage of some notable work on high-throughput surface characterization performed in the context of the Materials Project, for which the Pymatgen code has generated an efficient and user-friendly platform for generating surface slabs of inorganic crystalline materials.[40,41] Integrating this platform with open source software for graph theory applications, we apply a technique known as the max-flow min-cut algorithm[6] coupled with an advanced recursive implementation[42] to calculate the minimum cut for any particular surface. Examining the statistics of this minimum cut across all IZA zeolites yields a simple structural descriptor that predicts whether the material is likely to be synthesized as a layered 2D zeolite. Furthermore, the formulation of this graph theory problem is flexible and can be used to mimic the chemistry of specific 2D zeolite synthesis methods such as the ADOR strategy[23] by reweighting edges in a graph to model varying bond strengths. Thus, we show that we can mathematically formalize the generation of zeolite surfaces using graph theory and use this information to make useful predictions regarding the synthesizability of 2D zeolites.

Methods: Theory of Minimal Cuts in Graphs

Motivation

We are interested in the fundamental question of how a 3D material can be cut into two separate partitions by removing a minimum number of bonds or, given pairwise potentials, the minimum total bond energy per unit area. Our assumption is that this minimum cut indicates that the surface formed is energetically preferred[43] or that the delamination of the 3D structures into 2D sheets may be facile in this direction. While the true surface energy under synthesis conditions is extremely complex due to solvent, pH, structure directing agents (SDAs), etc., this minimum cut solution serves as a simple starting point from which an interesting structural descriptor will later be derived. However, we note that conceptual use of minimum cut surfaces has been able to predict and rationalize the crystal dissolution of a zeolite by atomic force microscopy,[4] and more recently it has been used in coarse grain crystal growth modeling.[5] Before entering a detailed discussion of the methods, the solution to this minimum cut problem for zeolites is visualized in Figure for zeolite MWW where the minimum cut surface termination of the (001) and (100) surfaces is determined. For example, the (100) surface cannot be formed by removing fewer than two bonds per unit cell.

Figure 1

Zeolite MWW is shown along with the min-cut surface terminations for both (001) = blue and (100) = green Miller planes. For either Miller plane, the surface termination shown represents the minimum number of bonds that can be cut while preserving the 2D periodicity corresponding to that particular Miller surface. Color code: yellow – silicon, red – oxygen, white − hydrogen.

Within graph theory, the problem of determining the minimal cuts (see next section for formal definition) arises in many varieties and has been an active area of research for decades[44−50] with importance in a large number of applications.[51−55] We will combine previous graph theory work on minimal cuts[42] in a novel application: the automatic determination of minimum or near-minimum cut surface terminations for any given Miller face of any given zeolite. Since a zeolite can be mathematically interpreted as a simple connected graph, we can use a minimum cut algorithm to solve the surface termination of a particular Miller face that cuts the minimal number of Si–O bonds. In this Methods section we start with mathematical notation and briefly outline the minimal amount of information necessary to understand a special type of graph cut, known as the minimum s-t edge cut, and provide references for additional details. After explaining how a “zeolite graph” can be created, we show how the minimum s-t cut problem can be applied to minimize the number of bond cuts necessary to generate a particular Miller surface of a zeolite.

Preliminaries

The minimal amount of formalism is presented here in order to define the graph theory concepts utilized in this work, and the exact notation of ref (42) is used. We take G = (V, E) to be an edge-weighted, connected, directed graph with a set of vertices (or nodes) V and a set of edges E ⊆ V × V where each pair of vertices is ordered, and n = |V| is the order of the graph (number of nodes) and m = |E| is the size of the graph (number of edges). An edge is denoted by its vertex pair e = (u, v). The weight of an edge e = (u, v), denoted w, is a numerical value associated with that edge, and the weights of all edges are given by the list w = (w,w, ..., w). In this work the only possible edge weights will be w = 0, w = 1, or w = ∞, as explained later. Two different vertices in V may be specially distinguished as the source and target vertices, or s and t, respectively.

Now, a directed s-t path in G is any path which starts at s and ends at t, or more formally a sequence of vertices and edges of the form s, (s, v1), v1, (v1, v2), v2, ..., v(v, t), t. A generic s-t edge cut is a set of edges C belonging to E that, when removed from the graph, interrupts all paths from s to t. The value or weight of the cut, w(C) = ∑w, is simply the sum over the weights of all edges in the cut. From here on we will use the definition of a minimum cut (or min cut), denoted as C0, to describe an s-t edge cut whose weight w(C0) is a minimum among all possible s-t edge cuts. Henceforth we may drop the s-t for brevity since all graphs in this work have a source and target. Multiple minimum cuts can exist, in which case C0 is a set of min cuts. Finally, a near-minimum cut (or near-min cut) Cϵ is an s-t edge cut whose weight is w(Cϵ) ≤ (1 + ϵ)w(C0), in which ϵ is a threshold to control how “near” a near-min cut must be to the min cut. Again, Cϵ is a set of near-min cuts if more than one exists. The following section gives visual examples of this formalism and briefly describes how the min cuts are calculated.

Determining the Min and Near-Min Cuts

The min cut of a directed, single-source, single-target graph is determined by computing the maximum flow of the graph.[56] Several algorithms can identify the max flow, and in our work this problem is solved using the “shortest-augmenting path” as implemented in the python library NetworkX.[57] Once the max flow is known, the value of the min cut is known by the max-flow min-cut theorem, and a single C0 solution can be easily determined.[6] Further details can be found in the extensive literature regarding the max-flow min-cut theorem, and we point the reader to ref (42) and references therein for more information. We show an edge-weighted, connected, directed graph with source s and target t in Figure to illustrate the min (C0) and near-min cuts (Cϵ) that can be computed. Note that each black double-sided arrow represents two antiparallel edges which effectively remove the directionality of all s-t paths in the graph. This is desired because materials’ bonds have no directionality but the min-cut algorithm in this work operates only on directed graphs. As a consequence the min-cut value is modified to be w(C) = ∑w/2 when antiparallel edges exist to avoid double counting.

Figure 2

Sample graph with source and target nodes (purple) for which we seek to find all min and near-min cuts, where red edges represent those included in the cut. Panel a shows the min-cut solution while panels b and c show near-min-cut solutions. Hence C0 = (C) and Cϵ = (C, C, C) when ϵ = 0.125. More near-min cuts could be included in Cϵ with an increased value of ϵ.

Determining the first min cut (Figure a) is relatively easy using the max-flow min-cut theorem. However, identifying all min and near-min cuts (Figure a–c) requires additional effort, for which Balcioglu and Wood have proposed an elegant algorithm that can be easily implemented.[42] First, a single min cut C0 is found by the max-flow min-cut theorem. Next, a recursive call to the max-flow min-cut algorithm is performed, but at each level of the recursion tree, the weight of a particular edge is modified to force its inclusion or exclusion from the previously identified min cut. We refer the reader to Section 2.2 and Figure 2.2 of ref (42) for more specifics and to the Supporting Information for the Python implementation of this recursive function. This recursive search of possible cuts ultimately outputs all min cuts and any near-min cuts based on the user specified value of ϵ. It should be noted that the computational feasibility of identifying Cϵ decreases for increasing ϵ since the number of cuts may be exponential in the size of the graph; however, we are only interested in values of ϵ = 0 in this work but still envision the use of nonzero ϵ in future work.

Creating an Initial Zeolite Nanosheet

A “naive” zeolite surface slab (nanosheet) can easily be generated using open source materials science libraries, and for this work we utilize the generalized surface slab building feature of Pymatgen.[40,41] The algorithm can build a surface slab of any Miller index from a bulk unit cell of any Bravais lattice, while advanced features can be used to build slabs where the two slab surfaces share an inversion point symmetry (when possible), to work with polar surfaces, and to determine symmetrically unique Miller faces. To generate our library of zeolite surfaces, we adopt the following procedure. We create a surface slab for a given Miller face of a given zeolite using Pymatgen and ensure that the slab contains Laue symmetry (when possible) so that both surfaces of the slab are symmetrically equivalent. A schematic representation is shown in Figure where, by Pymatgen convention, the slab’s surface is parallel to the ab plane of the new unit cell and the c-direction corresponds to the vacuum.

Figure 3

(Top) Each Miller surface slab contains N layers, where the slab thickness of the 2 ... N–2 bulk layers was chosen to be greater than 30 Å. (Bottom) The (103) Miller surface slab of MFI for N = 6 is shown with oxygen atoms excluded where the color of each Si corresponds to its layer in the slab.

We also require that the Pymatgen generated initial slab have N layers with translational symmetry where the thickness of the 2 ... N–2 layers must be greater than a certain cutoff distance. We chose a cutoff distance of 30 Å to minimize the self-interactions of the two surfaces so that the slabs could be used for first-principles calculations in future work. When enumerating all near-min cuts, we will only finalize a cut C0 or Cϵ if at least one Si atom attached to any of the edges e ∈ C0,Cϵ is in Slab Bulk Layer 1, shown in Figure . Imposing such a constraint eliminates the possibility of identifying equivalent cuts in the translationally symmetric layers, since a cut that only contains edges with Si atoms in Slab Bulk Layer 0 would be discarded due to the existence of an equivalent cut in Slab Bulk Layer 1. This also preserves the bulk character of Slab Bulk Layers 2 ... N–2 so that the surfaces remain sufficiently separated.

Application of Min Cuts to Zeolites

A zeolite slab can be easily converted into a simple graph upon which the min cuts can be calculated. O atoms are ignored since all are two coordinated to exactly two different Si. Thus, a zeolite can be interpreted as a simple connected graph where each node represents an Si atom with an edge between any two Si atoms that are connected by the same O. First, a zeolite slab is prepared and then interpreted as a simple graph, which is schematically shown in step 1 of Figure . The real zeolite graph should have edges that represent bonds crossing periodic boundaries in two dimensions (as denoted by the dashed lines in this toy example) but none in the third dimension, which corresponds to the direction perpendicular to the vacuum on either side of the slab. Any nodes in step 1 that are less than 4-coordinated are undercoordinated Si atoms whose bonds have been removed in the initial slab generation. They are identified as the initial surface nodes in step 2 and colored purple. On each side of the slab, these initial surface nodes are connected to a new single surface node, which is designated as the source or target node for a min s-t cut computation. The weights of these new edges must be set to infinity as shown in step 3 to ensure that any C0 and Cϵ are independent of the initial surface given in step 1. In step 4, we show a solution to the min s-t cut that is identified by the algorithm, specifically the leftmost cut. There is no need to repeat the algorithm to identify the analogous cut on the right side of the slab since it is immediately determined by Laue symmetry or, in its absence, translational symmetry between the slab layers which was visualized in Figure . The final step is to remove the partitions of the graph containing the source and target nodes, as well as passivate each dangling Si–O bond with an H atom. Hence we have identified a surface termination of the zeolite by disconnecting a minimum number of Si–O bonds since w = 1. It should be noted that the generated surface may express Si–OH, Si–(OH)2, or Si–(OH)3 groups to minimize the total number of removed bonds.

Figure 4

(Step 1) Schematic representation of a toy zeolite slab graph where nodes correspond to Si atoms and O atoms are replaced by a single edge. Note the dashed edges which are periodic when embedded in the 2D unit box. (Step 2) Undercoordinated surface nodes are identified. (Step 3) s and t nodes are connected to the surface nodes by edges with infinite weight. (Step 4) The min-cut solution is found. (Step 5) The dangling bonds in the final structure are passivated with hydrogen.

We also note that the algorithm for finding min cuts using the max-flow min-cut theorem operates on directed graphs, but in this section the zeolite graph in Figure is presented as a simple undirected graph. However, finding a min cut in undirected graphs is straightforward provided the following standard transformation is performed:[50,58] every edge in the undirected graph can be replaced with two directed antiparallel edges (Figure ), each with a weight equal to that of the original undirected edge. The previously discussed min-cut algorithm can be executed on this analogous directed graph, and any min cut represents a cut of equal weight in the original undirected graph, provided the cut weight is modified to w(C) = ∑w/2 to avoid double counting the additional antiparallel edge.

The Lammps Interface code,[59] which interprets nanoporous materials as periodic graphs using Python’s NetworkX package,[57] was used to convert zeolite structures into their analogous graphs, and Pymatgen was used to generate the initial surface slabs. We extended NetworkX’s default max-flow min-cut computation with the previously described recursive algorithm of ref (42) to solve the min cut of the initial slab. Using GNU parallel[60] to run embarassingly parallel jobs was sufficient to generate just one min-cut slab for each symmetrically unique Miller surface up to maximum index of 2 of every IZA zeolite in just a few days on an 8 processor desktop. This yielded a total of ∼3700 slabs.

Results and Discussion

Generation and Characterization of IZA Zeolite Surfaces

The min cut for each Miller surface up to a maximum index 2 for each IZA zeolite (excluding any interrupted structures, denoted with a “-” by the IZA commission) was solved to generate a library of 2D nanosheets. Thus, for each Miller face of each IZA zeolite we have calculated the w(C0) value, which is exactly equal to the number of edges that are cut to form the surface since w = 1 in this scheme. Since each slab is embedded in a new unit cell with the vacuum parallel to the c-direction as visualized in Figure , w(C0) is converted to a cut density by division with the area of the face spanned by the ab unit cell vectors. Then each Miller surface is ranked by this surface density of cut edges, δ = w(C0)/(|a × b|) with units of Å–2, as shown in Table for the examples of EMT and MWW. Here R1 and R2 denote the Miller surfaces with the lowest and second lowest min-cut densities, respectively. From now on δ will be referred to as the min-cut density, which can also be interpreted as a density of cut bonds when the weight of each edge in the graph is unity. Note that several Miller surfaces can have the same min-cut density due to symmetry equivalence, and this was exploited to save computational time by only performing the min-cut analysis on one of the symmetrically equivalent Miller planes as calculated by Pymatgen. Repeating the analysis in Table would be extremely arduous if not impossible by manual/visual inspection of all IZA zeolites. Some additional analyses on the statistics of all IZA min cuts are shown in the Supporting Information to highlight the necessity of the automated and robust algorithmic approach provided in this work.

Table 1

Ranking of the EMT and MWW Miller Surfaces Based on Their Min-Cut Density, δ [=] Å–2a

	EMT		MWW
rank	face	δ	face	δ
R1	(001)	0.0234	(001)	0.0112
R2	(100)	0.0248	(102)	0.0314
R3	(11̅0)	0.0248	(102̅)	0.0314
R4	(010)	0.0248	(11̅2)	0.0314
R5	(101)	0.0256	(1̅12)	0.0314
R6	(101̅)	0.0256	(012)	0.0314
R7	(11̅1)	0.0256	(012̅)	0.0314
R8	(1̅11)	0.0256	(100)	0.0331
R9	(011)	0.0256	(11̅0)	0.0331
R10	(011̅)	0.0256	(010)	0.0331

Surfaces ranked higher than 10 are omitted for clarity.

Application: Predicting IZA Zeolites Likely To Grow in Layered 2D form

We aim to have a predictor for zeolites which can grow in a stable, layered 2D form. A closer look at Table reveals a major difference between the statistics of the min-cut densities for MWW and EMT. For MWW, the difference in the min-cut density between the R2 and R1 surfaces, δR2 – δR1, is relatively large with a value of 0.02 whereas for EMT this quantity is practically zero. Now the question becomes whether such an outlying min-cut density of the R1 surface indicates that it can be more easily isolated during crystal growth, leading to the formation of layered 2D zeolites. In other words, can one more readily find synthesis conditions/SDAs to achieve enhanced stability of the R1 surface relative to other surfaces or obstruct growth in the dimension orthogonal to R1? To investigate this question, we plot δR1 vs δR2 – δR1 for all IZA zeolites to elucidate an important trend shown in Figure .

Figure 5

Plot of δR1 vs δR2 – δR1 where each data point corresponds to an IZA zeolite. Large values of δR2 – δR1 lead to a much higher probability that an IZA zeolite has a known 2D form, indicating it is a probabilistic descriptor for the ability to form such structures. The 15 IZAs with known 2D form and the top 15 IZAs with no known 2D form are listed in order of decreasing descriptor value to help identify them in the figure. The chemical formula given by the IZA Commission[35] (omitting counterions and SDAs) is also provided.

In this plot each data point corresponds to an IZA zeolite, where gold points have a known 2D form[22,25,61] and blue points have no known 2D forms. Picking a structure at random, one would have an approximate 7% chance of choosing an IZA zeolite that is known to exhibit a 2D form. However, when focusing on structures with the largest δR2 – δR1 values, that probability increases significantly (up to ≈70%). In other words, δR2 – δR1 is a probabilistic descriptor for identifying an IZA zeolite with a known 2D form. Physically, this simple structural descriptor says the following: it is more likely to find synthesis conditions or SDAs that block the growth of a zeolite in one crystallographic dimension when the face perpendicular to that dimension has a much lower minimum cut density than any other face. This leads to the formation of a 2D layered precursor, and it is only the growth in the third crystallographic dimension that is blocked. To ensure that the descriptor actually predicts a 2D zeolite with the correct surface, we manually investigated all reports of the known 2D structures from the literature following the references in refs (22) and (25). In all cases, the reported 2D structure[62−70] is formed such that the expressed surface corresponds to the same surface we identify as R1 in our high-throughput screening. To summarize, Figure clearly demonstrates that the R1 surface is much more likely to be expressed under given synthetic conditions than the R2 surface when δR2 – δR1 is large, while the absolute value of δR1 is less relevant for predicting layered 2D crystal growth.

We note that real zeolite synthesis is a very complex process, which requires fine control of reaction conditions, including but not limited to the right reactants, a specific structure directing agent, controlled pH (of reaction medium), optimized reaction temperature and time, and the right mixing fraction and procedure. The experimental realization of our predicted 2D zeolite candidates also relies on the fact that the right reaction conditions need to be identified. This is beyond the scope of our current work. However, our predicted surface terminations of the Miller planes with the lowest cut densities may provide crucial insights on one of the most important reaction ingredients of zeolite synthesis, i.e., the structure directing agents. By computationally screening the binding strengths of conventional and unconventional SDA molecules with different Miller planes of a selected 2D zeolite candidate, one may identify a SDA molecule with much stronger interaction with the lowest cut density Miller plane than with other Miller planes, and therefore zeolite growth will be promoted along this unique direction and inhibited along other directions, which results in a 2D zeolite with a maximally expressed Miller surface. We also note that the lifetime of zeolite surfaces has been reported to inversely correlate with the density of surface dangling bonds (i.e., the cut density),[4] indicating that the lowest density cut Miller planes will have more time to interact with SDA molecules and have increased expression in the growth process, furthering the chance that a 2D instead of a 3D zeolite will be formed. It should be noted that MFI represents an outlier since it is the only zeolite that can be synthesized in a stand-alone 2D form,[71] i.e., not as a layered 2D precursor, and also has δR2 – δR1 approximately equal to zero. This 2D form is also achieved by sythesis methods unique to this structure, namely, as a multilamellar precursor with surfactant.[22] Thus, a large δR2 – δR1 value more reflects the probability to find synthesis conditions to direct formation of 2D layered precursors but as expected does not represent the complex surfactant chemistry at the MFI surface that is utilized to direct its stand-alone 2D growth.

Finally, the structure corresponding to each data point is provided in the Supporting Information. These results provide a list of structures that can be immediately targeted experimentally because they are higher probability candidates to form novel 2D layered structures according to the statistics of currently known 2D layered zeolites. These large δR2 – δR1 value structures serve as a starting point for future experimental and computational efforts to predict and investigate which synthetic conditions and SDAs may result in some 2D layered structure where the expressed surface is defined by our predicted R1 surface.

Application: Potential 2D Zeolites for Water Desalination

The nonequilibrium molecular dynamics simulations of Jamali et al. demonstrated that using 2D zeolite nanosheets for water desalination could hypothetically provide large improvements in water permeation performance over current technologies.[31] We now look at all known zeolites with 2D form (as well as those in the “high probability zone” of Figure ) to determine which structures have potential for this application, i.e., are porous in the direction perpendicular to the R1 surface. Here we define a nanosheet as porous if the largest free sphere (or pore limiting diameter), Df, in the crystallographic direction perpendicular to the R1 surface is larger than the kinetic diameter of water. Df in the direction perpendicular to the R1 surface was calculated with Zeo++[18] and is plotted vs the descriptor in Figure .

Figure 6

Largest free sphere, Df, in the direction perpendicular to the R1 surface plotted vs the descriptor δR2 – δR1. For separations to occur in such 2D structures, this value must be larger than the kinetic diameter of the smallest species in the separation mixture. MFI is the only known 2D zeolite to achieve porosity through the nanosheet, but nine potential structures are highlighted that would also be porous and have not yet been discovered in a 2D form. The chemical formula given by the IZA Commission[35] (omitting counterions and SDAs) is also provided.

Only MFI satisfies this porosity criterion for all zeolites with a known 2D form. SOD and MWW have the next largest values for Df, but neither has pore limiting diameter large enough in the direction perpendicular to the R1 surface to allow water to pass from one side of the nanosheet to the other. However, several structures can be identified with relatively large δR2 – δR1 value (indicating a higher probability of having a layered 2D form) that have a pore limiting diameter through the R1 surface larger than the kinetic diameter of water, some of which are listed in Figure . Our analysis shows that these materials have structural characteristic similar to those of other stable 2D layered zeolites, and, given their potential for different types of separations, we hope that this work will stimulate a more targeted effort to synthesize them.

Application: Predicting Likely 2D Zeolites from ADOR Disassembly

There exist a variety of specific synthesis techniques for achieving a 2D form of an IZA zeolite.[22,25] However, if one were to ask which materials may form layered 2D sheets for a specific synthesis method, it is possible in some situations to tailor the min-cut algorithm to mimic the chemistry of a particular synthesis technique. Here we provide an example on the versatility of our algorithm to show how we can identify potential 2D zeolites formed during the disassembly step of a special synthesis procedure, namely, the ADOR strategy. In this technique, germanium preferentially occupies double 4 ring (D4R) sites which are selectively hydrolzyed upon acid treatment. With this knowledge, the weights of the edges in a zeolite graph should be set to zero if the edge contains one node in a D4R unit and one node outside the D4R unit. Removing the penalty to cut these edges mimics the selectivity of O–Ge bond hydrolyzation. This reweighting of bonds attached to D4R units is shown schematically in Figure a.

Figure 7

(a) A schematic showing the reweighting of edges corresponding to bonds that are selectively hydrolyzed during ADOR disassembly. (b) The Df vs our descriptor, δR2 – δR1, is plotted for IZA zeolites containing D4R motifs. Only structures where δR1 = 0 are considered, so a nonzero value of the descriptor indicates that exactly one Miller surface (the R1 surface) can be generated by only cleaving bonds connecting a D4R unit to the rest of the structure. The red line once again indicates the kinetic diameter of H2O.

With this new weighting scheme, a min-cut density of zero, or w(C0) = 0, can occur for a particular Miller plane if all bonds in the min cut correspond to those that are selectively hydrolyzed during ADOR disassembly. Thus, the first requirement for 2D sheet formation with this technique is that δR1 = 0. However, if there exist two min-cut surfaces with δ = 0 for two different Miller planes, it is evident that no 2D sheet could form as it would be hydrolyzed into discrete fragments that lack 2D periodicity. For 2D layers to be formed, it is sufficient to see that one and only one Miller surface has a min-cut density of zero, or in other words, the second requirement is that δR2 – δR1 ≠ 0. Figure b shows the pore limiting diameter vs δR2 – δR1 when the weighting modification of Figure a has been applied to demonstrate promising materials for ADOR disassembly.

Here only IZA structures that have D4R units are shown and are color-coded by those that have been synthesized as germanosilicates. Reference (23) determined, presumably by manual inspection, that “the most suitable candidates for top-down synthesis of 2D zeolites [are]: ITG, ITH, ITR, IWR, IWW, SVV, UOS, and UTL ... [and] IWV”, and our automated approach provides very similar prospective. Regarding nongermanosilicates, ref (23) proposes IWV as a potential material for ADOR disassembly if it can be synthesized as a germanosilicate, but excludes IFY, UFI, UOV, and ISV (which were likely excluded since they contain D3R or single 4 rings that, if hydrolyzed, would destroy the 2 dimensionality). However, we additionally highlight as a target for the synthesis community that DFO can be included as a potential candidate if the germanosilicate version of the framework can be synthesized. UOS appears in ref (23) but not our list since our calculations revealed that more than one Miller surface can be formed by only hydrolyzing bonds connected to D4R subunits, or δR2 = 0. Clearly ITG, ITR, IWR, and IWW are all high potential candidate structures for separations applications since the 3D → 2D transformation via Ge–O hydrolysis results in an atomically thin sheet porous to water molecules. Finally, it should be highlighted that this reweighting scheme can generally be applied to any structural motif in zeolites (e.g., double 6 rings, cages, etc.) to accommodate future disassembly techniques. Subsequent min-cut calculations can then provide a formal way of determining whether a 2D sheet can be formed from a 3D structure by only breaking bonds connected to specific building blocks.

Conclusions

We have applied a powerful graph theory technique to find and enumerate the minimum cut 2D surfaces of 3D crystallographic materials. The ability to calculate a Miller surface termination that minimizes the number of cut bonds (or the minimum total energy of cut bonds given a classical potential) provides a mathematical approach to identifying important surfaces, and we have specifically applied these ideas to zeolites in this work. For any given zeolite and Miller face, one can automatically and rigorously enumerate all of the minimal and near-minimal cuts to create a library of 2D nanosheets. To our knowledge this is the first in silico design approach using graph theory to study crystal surfaces in such a high-throughput manner. While specifically applied to IZA zeolites in this project, this methodology has the potential to be applied to other crystallographic systems (i.e., ZIFs, metal oxides, MOFs, etc.)[72−75] to investigate various surface terminations in a formalized, high-throughput methodology. One could also, for example, perform an identical analysis for the database of aluminophosphates[76] and other zeotypes due to the generality of our graph theory approach.

Our algorithmic approach to zeolite surface generation not only yields a probabilistic descriptor for the likelihood that an IZA material has a known layered 2D form but also correctly identifies the expressed surface. From random selection one has ∼7% chance of selecting a material with known 2D form. Using our descriptor, one can bias this selection probability to ∼70%. This indicates that materials with favorable descriptor values and as of yet unknown 2D form are the most likely to be discovered in layered 2D form upon investigation of new synthesis conditions. We provide a list of structures which can guide experimental efforts for attempted synthesis of layered 2D zeolites with far higher probability than random search. We have furthermore demonstrated the versatility of the algorithm by predicting suitable candidates for a particular synthesis method, namely, the 3D to 2D transformation during the disassembly step of ADOR. While ADOR is relatively new, Eliášová et al. commented that new methods to selectively design structural weakness at specific structural motifs (other than Ge–O in ADOR) and to selectively break such bonds will be critical to developing new 2D forms of known zeolites.[23] When such synthetic procedures are discovered and developed, our algorithmic approach will be invaluable to automatically identify likely materials for such 3D to 2D transformations from the large number of IZAs, non-IZAs, and hypothetical zeolites. Assuming sufficiently accurate classical potentials, a natural extension of this combined ADOR/min-cut analysis would be to identify other inorganic crystals that are high probability candidates for exfoliation, or in other words exhibit a 3D to 2D disassembly with exactly one preferred Miller surface.[75]

Our formal approach to enumerating zeolite surfaces also opens a path for other computational studies that can be performed to better investigate and understand zeolite surfaces. For example, screening based approaches can be applied to identify SDAs or solvent conditions that energetically favor the min-cut surface over higher density cut surfaces, leading to controlled structure growth. Informatics based identification of SDAs may be able to predict compounds that will lead to isolation of the 2D layered form of some of the high-potential candidates identified in this work. Ultimately, the theory applied here will not only be important for the continued investigation of zeolite surfaces and identification of potential 2D zeolites but also will provide a new methodology to examine surfaces of other classes of crystallographic materials.