## Abstract

Zeolites are microporous crystalline aluminosilicate materials whose atomic structures can be usefully modelled in purely mechanical terms as stress-free periodic trusses constructed from rigid corner-connected SiO_{4} and AlO_{4} tetrahedra. When modelled this way, *all* of the known synthesized zeolite frameworks exhibit a range of densities, known as the flexibility window, over which they satisfy the framework mechanical constraints. Within the flexibility window internal stresses are accommodated by force-free coordinated rotations of the tetrahedra about their apices (oxygen atoms). We use rigidity theory to explore the folding mechanisms within the flexibility window, and derive an expression for the configurational entropic density throughout the flexibility window. By comparison with the structures of pure silica zeolite materials, we conclude that configurational entropy associated with the flexibility modes is not a dominant thermodynamic term in most bulk zeolite crystals. Nevertheless, the presence of a flexibility window in an idealized hypothetical tetrahedral framework may be thermodynamically important at the nucleation stage of zeolite formation, suggesting that flexibility is a strong indicator that the topology is realizable as a zeolite. Only a small fraction of the vast number of hypothetical zeolites that are known exhibit flexibility. The absence of a flexibility window may explain why so few hypothetical frameworks are realized in nature.

## 1. Introduction

In the mid-eighteenth century, Axel Fredrick Cronstedt discovered that when certain minerals were heated in a flame they ‘emit gas and puff up’ [1,2]. Accordingly, he assigned the name *zeolite* by merging the Greek words *ζειν* and *λιθoσ*, meaning ‘boiling stone’. In his short article, Cronstedt made the prescient remark that a ‘more copious supply of pure material, which is not presently available, would enable us to ascertain whether these substances are of any utility in any handicraft’.

It took over two centuries before the remarkable utility of zeolites could be realized. It first had to be understood that zeolites are microporous, containing periodic arrays of tunnels and pores in two and three dimensions that are large enough to allow small molecules, such as water and light hydrocarbons, to diffuse relatively unimpeded. Second, it had to be learned how to synthesize a ‘more copious supply of pure material’ before they could be exploited with commercial success [3–5]. The theoretical occupiable internal surface areal density can exceed 2000 m^{2} cm^{−3} in some zeolite structures. Further, the chemical composition of zeolites can make them behave as solid acids, which is a valuable property for cracking hydrocarbons. Their ability to behave as molecular sieves and catalysts, and remain stable at high temperatures, makes zeolites an important component of the modern-day petrochemical and fine chemical industries.

There are presently 206 distinct zeolite framework types known [6,7]. About 50 of these occur as minerals, and 150 have been found through exhaustive synthesis efforts over the past 60 years, and have yet to be found in nature. Most new framework types are discovered by hydrothermal synthesis, and each discovery of a new framework topology is essentially serendipitous.

There has been a parallel effort to discover new zeolite framework topologies by mathematical modelling. A number of ingenious strategies have been devised for generating new frameworks, many of which have been summarized in [8]. The primary goal of the modelling has been to discover new framework types whose computed properties are of interest to the petrochemical industry. A secondary goal relates to the fact that these structures can be fascinating, and even beautiful, and we can learn much about crystal framework topology by studying them. There now exist databases containing many millions of hypothetical frameworks [9,10], but we still do not know how to guide the synthesis of the many ‘designer zeolites’ found in these databases.

As the size of hypothetical zeolite databases has expanded, it has become a puzzle as to why there are so many hypothetical candidates, yet so few realized materials. This large discrepancy between what mathematics proposes and what nature provides has been referred to as the ‘zeolite conundrum’ by Blatov *et al.* [11]. Although the number of hypothetical structures is theoretically infinite, it seems unreasonable to assert that therefore there is an infinite number of realizable zeolites. It is likely that the databases have been overly inclusive. Our database [12], for example, is built on the assumption that the materials comprise pure silica, SiO_{2}, and the framework energies are computed under this assumption. In order to be inclusive of other compositions, a generous spread in framework energy is allowed. This, undoubtedly, allows many unrealizable frameworks, as silicates, to be included in the database. A feasibility factor has been proposed for pure silica frameworks [13] based on a computation of the framework energy of formation using the Sanders–Leslie–Catlow potential in the GULP program [14]. It has also been observed that about 80% of realizable zeolites can be decomposed into ‘natural building units’, which are simple three-dimensional tiles or cages. However, there are too many exceptions to these rules for them to act as decisive framework discriminators. The question is, ‘Are there other topological or physical properties that determine framework feasibility?’ We believe that framework flexibility is that key property and we explore it here.

Although zeolites are crystalline assemblies of atoms, they can also be usefully viewed as framework oxide materials that are closely related to quartz. An excellent working approximation is that a zeolite is a tetrahedral silicate, of average composition SiO_{2}. At the local level, the structure can be thought of as a periodic array of corner-sharing SiO_{4} tetrahedra, with shared oxygen atoms being located at the tetrahedral apices, and silicon atoms nestled inside the tetrahedron (figure 1). Computer modelling, with realistic atomic force fields, shows us that compressive strain tends to be absorbed mostly by coordinated rotations of the tetrahedra about the oxygen linkages. Conversely, large tensile strains tend to be absorbed by deformations of the tetrahedra, particularly stretching of the Si–O bond distance, and changes of the O–Si–O angle away from the ideal tetrahedral angle, .

Zeolites can also be modelled by an idealization, wherein the tetrahedra are considered to be perfect and rigid, and the corner linkages (oxygen atoms) act as force-free spherical joints. Although this idealized force field represents a radical simplification of the subtle crystal chemistry, it allows us to probe some of the fundamental properties of zeolite structures. For example, Dove and co-workers [15–17] have used such a model to explore rigid-unit modes (RUMs) in tetrahedral frameworks and their role in phase transitions in these materials. As we discuss in more detail here, it also allows some of the basic properties of zeolite flexibility to be explored.

It has long been known that some zeolite frameworks are flexible, sodalite [18] and zeolite rho being prominent examples. However, it was realized only recently that *all* of the known zeolite frameworks are flexible and exhibit a range of densities over which they are stress-free when represented as ideal frameworks [19–21]. This range of densities is referred to as the *flexibility window*, and the flexibility index is defined as the fraction [20,21]. A preliminary table listing the flexibility indices of most zeolites was presented in [21]. Our methods for finding flexibility modes have improved since that report, and the handful of frameworks that were reported as lacking a flexibility window, most notably the mineral goosecreekite, have since been found to be flexible after all when represented in lowest symmetry, *P*1, and when allowing mixed Si and Al tetrahedra.

The mathematical signatures of flexibility in a framework come from rigidity theory, which is concerned with finding necessary and sufficient conditions for generic rigidity of frameworks. While the question of extending rigidity theory to the crystallographic setting had been open for several decades, major breakthroughs have been made by Malestein & Theran [22]. In the ideal model of a zeolite, the flexibility window represents a zero-energy nullspace of possible framework configurations. Curiously, neighbouring and overlapping points in the flexibility window are frequently remote from each other topologically; unrelated folding mechanisms can have identical densities. In general, the only sure way to transform one folded conformation into another is to expand the framework to its minimum density, and then refold the framework via a different mechanism. However, there are no energy gradients to guide the fold. Consequently, our computational tools that attempt to demarcate the full extent of the flexibility window have relied on stochastic methods rather than gradient-guided exploration.

In this paper, we describe some of the progress that we have made in exploring mechanisms within the flexibility window. The computational tools we use are based on the singular value decomposition (SVD) methods developed by Pellegrino and co-workers [23,24], Guest & Hutchinson [25] and Hutchison & Fleck [26] for exploring mechanisms allowed by the compatibility matrix of a three-dimensional periodic truss. Essentially, we are treating zeolites as ‘deployable’ periodic frameworks. We express the compatibility relations in dynamical form, and derive an expression for the ‘vibrational’ entropy, that is valid at high temperatures. The null modes have zero vibrational frequency, and thus potentially dominate the entropy, contributing a configurational entropy component. In a previous paper, we explored the configurational entropy in frameworks allowing mechanisms to adopt low symmetry [27]. Here, we explore the configurational entropy while maintaining high symmetry and find that the configurational entropy can also possess a strong component from low-frequency soft modes.

## 2. Zeolites as three-dimensional periodic trusses

We approximate the zeolite framework by idealizing it as a three-dimensional periodic network of corner-connected rigid tetrahedra. This is accomplished by imposing a system of distance constraints between oxygen atoms (referred to as ‘O-atoms’, each of mass *m*) that define rigid tetrahedra (figure 2). For a unit cell containing *N* O-atoms, there will be 3*N* O–O distance constraints. This is because each O-atom is a vertex in a degree six graph, that is, there are six bars connected to each oxygen joint. Each edge is shared by two vertices (O-atoms). Silicon atoms (referred to as ‘T-atoms’) sit at the geometric centre of each tetrahedron. For simplicity, we ignore the T-atoms because they only contribute a passenger mass to the truss, which is now defined by the sixfold coordinated O-atoms. This O-atom-based description of a zeolite framework is not the one that is usually used when describing zeolite structures. It is more common for a zeolite framework to be described as a degree four graph with T-atoms located at the vertices and the O-atoms being ignored as passenger atoms near the midpoints of edges (figure 2).

In real silicate materials, it is the Si–O chemical bond that dominates, and essentially defines, the structure. The SiO_{4} (or TO_{4}) tetrahedra are usually found to be near perfect, with O–Si–O angles lying very close to the ideal tetrahedral angle of . The Si–O–Si angles in zeolites are less constrained and tend to range between 130° and 180°, with an experimental mean angle of 154° [28]. As an idealization, we treat the intertetrahedral T–O–T angle as being unconstrained, with each O-atom behaving as a force-free spherical joint. This idealization of the zeolite framework, as a periodic truss of rigid tetrahedra, has led to several insights into their structural properties. Zero-frequency RUMs have been discovered, and their role in phase transitions and the absorption properties of framework materials has been investigated [15,16,29–31]. It has also been discovered that all of the known zeolite frameworks are flexible, exhibiting a range of densities over which volume changes are absorbed by tetrahedral rotations, without generating any stresses [19,21]. The flexibility window represents the nullspace of the framework elastic energy that can be traversed by many unique folding modes. Rigidity theory provides us with powerful tools for exploring these folding modes.

From rigidity theory, we know that
2.1where **C** is the compatibility matrix, **d**_{v} is a column vector representing infinitesimal displacements of each vertex and **e** is a column vector representing the infinitesimal elongation of each joint. Here, we use the Cartesian coordinate system. Equation (2.1) is valid only for infinitesimal displacements because **C** itself evolves as the vertex coordinates change. In principle, for fixed unit cell dimensions of a system containing *N* tetrahedra, **C** can be represented as an *N*×*N* tensor whose elements operate on the *N*×1 column vector containing the individual infinitesimal vector displacements **d**_{i}=(*δx*_{i},*δy*_{i},*δz*_{i}). For *N* tetrahedra, there will be 2*N* shared joints. For computational purposes, it is easier to expand **d**_{v} and **e** as 3*N*×1 column vectors that contain *δx*_{i}, *δy*_{i} and *δz*_{i} as separate scalar components. This, in turn, requires **C** to be represented as a 3*N*×3*N* matrix.

In a rigid framework, the edge extension field is zero, **e**=0. In this limit, the elements of the compatibility matrix can be determined from the length constraint equation relating the distance between joint *i* and joint *j*,
2.2(*x*_{i},*y*_{i},*z*_{i}) are the Cartesian coordinates of joint *i*. The *X*_{ij}, etc. are the periodic Cartesian translations required when one of the connected *i*, *j* pairs is a translated image of a joint on the other side of the unit cell. The *X*_{ij},*Y* _{ij},*Z*_{ij} ensure the periodic continuity of the topology. There are *N* unique constraint equations, one for each bar. Differentiating, and cancelling the factor of 2, we get,
2.3This is equivalent to one row of the expanded matrix equation **C****d**_{v}=0, relating to the edge *i*−*j*. For every joint *i*, there are four connected joints *j*. Thus, each row of the compatibility matrix **C** has four non-zero entries of the type (*x*_{j}−*x*_{i}−*X*_{ij}) in each column for the connected joint *j*. Each *i*,*j* pair, of course, has three separate components for the *δx*_{i}, *δy*_{i} and *δz*_{i} parts, and so we develop a 3*N*×3*N* matrix equation.

A zeolite framework with a fixed unit cell is *locally isostatic*. This means that, ignoring any symmetry, the number of degrees of freedom available to the framework equals the number of constraints acting on the framework. This can be understood by examining each tetrahedron, which has three translational and three rotational degrees of freedom, to make six degrees of freedom per tetrahedron. However, the *x*,*y*,*z* coordinates of the tetrahedral apices, where the O-atom joints are located, are shared by two tetrahedra. This means that there are 12 shared constraints, giving 12/2=6 constraints per tetrahedron. Each tetrahedron has, in general, a ‘balance’ of six degrees of freedom and six constraints, hence the term *locally isostatic*. As pointed out by Guest & Hutchinson [25], if we allow the unit cell dimensions to vary, this adds six degrees of freedom to the periodic system. Three are associated with the cell edge lengths *a*, *b* and *c*, and the remaining three with the interaxial angles *α* (between ** b** and

**),**

*c**β*(between

**and**

*c***) and**

*a**γ*(between

**and**

*a***). Three of these degrees of freedom resolve themselves as the trivial linear translations of the whole lattice. However, the other three are not the trivial rotations, as these would render the periodic framework structure incommensurate with the underlying lattice. Instead, the remaining three resolve themselves as internal shear modes [25]. These three internal modes guarantee that the framework has at least three infinitesimally flexible modes [25]. The additional degrees of freedom introduced by variable unit cell dimensions can be important to the framework flexibility, and need to be included in the flexibility analysis.**

*b*The unit cell volume is given by
2.4If the unit cell volume is to be held constant, this reduces the additional degrees of freedom to five, and the interdependency between the six cell displacement parameters, *δa*, etc., can be found from equation (2.4) by asserting *δV* =0.

Once we allow the unit cell to change, the components of the compatibility matrix are no longer as simple as we have described above, because changes in the cell parameters and changes in the joint coordinates are not linearly independent. Any change to the cell parameters affects the Cartesian locations of all the joints, which in turn affects all of the 3*N* elements *δx*_{i}, etc. in the displacement vector **d**_{v}. Care must be taken to account for the possibility that opposite ends of a joint may be in different images of the periodic cell.

To accommodate displacements related to cell changes, we introduce a modified displacement vector **d**, which is a (*N*+6)×1 column vector. The first 3*N* elements are the elements of **d**_{v} (which are the Cartesian components of all the relative displacements, ). The last six elements are *δa*, *δb*, *δc*, *δα*, *δβ* and *δγ* cell change parameters. The two are related via a 3*N*×(3*N*+6) transformation matrix **T**
2.5where the elements of **T** are obtained, three rows at a time for each vertex, from equation (A 4) in appendix A.

The tensile force field **t** on the O–O linkages (bars) is given by
2.6**G** is a diagonal 3*N*×3*N* spring constant matrix,
2.7*k*_{l} is the spring constant for member *l*=*ij*≡*ji*. The off-diagonal terms are zero because, in our idealization, we consider only elastic forces along each rod, and here we ignore angular terms between different rods. These terms are readily included but are ignored here in the interests of simplicity. These stresses translate into a 3*N*×3*N* force field **f** on the O-atoms (joints)
2.8The equilibrium matrix is simply **A**=**C**^{T} through the static–kinematic relation, where **C**^{T} is the transpose of **C**. Substituting for **t**, we have
2.9We now invoke Newton's second law of motion to obtain
2.10or, in terms of the modified displacement vector **d** (equation (2.5)),
2.11which is the simple harmonic equation. The matrix **C**^{T}**GC** has dimensions 3*N*×3*N*, and is a measure of the framework stiffness. **M** is also a 3*N*×3*N* diagonal mass matrix, and we are assuming, for simplicity, that the silicon atoms are massless.

In the presence of external applied forces, which we decompose into a bar-tension field **t**_{A}, the equation of motion becomes
2.12**d**_{A} is a displacement amplitude associated with the externally applied tensions. Of course, this dynamical equation applies only to infinitesimal displacements because the compatibility matrix, **C**, coevolves with the solution.

In the absence of external applied forces (zero pressure), **t**_{A} and **d**_{A} are both zero. The spring matrix can also be expressed in simpler form, **G**=*k***I**, where *k* is the spring constant of all the O–O bars, which are all of length *d*_{OO}, where **I** is the 3*N*×3*N* identity matrix. We try standard solutions of the type to obtain the eigenvalue equation
2.13which is valid for infinitesimal displacements about the nullspace configuration.

We introduce a wavevector ** K**, which has components (1/(

*n*

_{x}

*a*),1/(

*n*

_{y}

*b*),1/(

*n*

_{z}

*c*)), where the

*n*

_{x},

*n*

_{y}and

*n*

_{z}are the repeat lengths along the

*a*,

*b*and

*c*axes for a wave that is tilted by

**. Bloch's theorem tells us that for a periodic system any periodic function**

*K**ψ*(

*x*) obeys 2.14In this way, we can explore the null mode states within the Brillouin zone, and not just at

*K*=0. To designate the

*K*-dependence, we annotate the compatibility matrix as

**C**

_{K}.

Our equations of motion provide a way of estimating the configurational entropy of the flexible framework. Thermodynamically, heat provides the energy that allows the system to explore its neighbouring states. In the limit of high temperature *T*, the vibrational entropy of a system of uncorrelated oscillators is
2.15*N*_{D} is the number of degrees of freedom, which is *N*+6 if the trivial translation vectors are to be included. As usual, the Boltzmann constant *k*_{B}=1.381×10^{−23} J K^{−1} and the reduced Planck constant J⋅s. The angular frequency *ω*_{l}=2*πν*_{l} rad s^{−1}, where *ν*_{l} is the vibrational frequency of mode *l*.

The product of the eigenvalues of any matrix equals its determinant. Thus,
2.16Taking the natural logarithm, the configurational entropy can be expressed as
2.17with *S*_{0}(*T*) being a temperature-dependent offset. A more generalized form, retaining the elastic and mass matrices, is
2.18

The expression that is most useful to us for the periodic zeolite frameworks is equation (2.17), which is dominated by the zero-frequency null modes. The zero eigenvalues need to be handled by a renormalization method. We assign a constant small number, 1×10^{−5}, to all eigenvalues. Thus, each null mode adds 5*k*_{B} per unit cell to the entropy. As the premise is that motions are infinitesimal, our expression for entropy will be inexact for the null modes where there is a divergence. Such divergences imply that near-infinite volumes are being explored, which is not the case. Qualitatively, the singularities flag those configurations that have high densities of null modes. Consequently, our expression for configurational entropy is useful for comparative purposes.

The solutions (equations 2.17 and 2.18) incorporate an elastic potential term, which rigidifies when we set the spring constant . This creates an additional renormalization issue as . This can be dealt with by just examining the configurational entropy relative to *S*_{0}.

Although the Bloch wave approach provides a powerful means for evaluating the density and symmetry of null modes within the Brillouin zone, it generally fails if we actually follow a mechanism by folding the framework according to that mode. The folding mode for the general *K* vector will generally involve large supercells, and will immediately violate the lattice periodicity that underpins the Bloch wave approach. To bypass this problem in our analysis, we adopted a two-step approach. First, we used the Bloch wave formalism to examine zeolite frameworks in their high-symmetry (i.e. low-density) configurations. Then, to follow mechanisms corresponding to wavevector *K*=(1/*n*_{x},1/*n*_{y},1/*n*_{z}), we create a supercell by repeating *n*_{x}, *n*_{y}, *n*_{z} cell lengths along the *x*, *y* and *z* axes. Transformed to this supercell, we then need to follow only the *K*=0 modes. Depending on the instantaneous symmetry of the framework allowed by the mechanism, the number of modes can grow as (*n*_{x}*n*_{y}*n*_{z})^{p/3}, where *p*=0,1,2,3. Frameworks for which *p*=3 are extensive, meaning that the number of mechanisms, and hence the thermodynamic configurational entropy, grows linearly with the supercell volume.

We have written a code, the Zeolite Null Space Explorer (ZeNuSpEx), that implements the above algorithm. The null modes are found using SVD as advocated by Pellegrino & Calladine [23]. Efficient SVD algorithms are available for sparse matrices that can run rapidly on graphics-processing units [32]. The (*N*+6)×(*N*+6) square matrix has the same null eigenvectors as **C**_{K}, and so the flexible null modes may also be found by conventional, square-matrix, eigenvalue methods.

Being degenerate, the null modes are not necessarily orthogonal, and are not constrained to preserve any internal crystallographic symmetry. Although we did not implement it, the eigenproblem could be simplified considerably by imposing symmetry constraints on the elements of the displacement vector **d**. The length of **d** and **T** can be reduced so that they act only on the unique (basis) elements of the displacement vector, thereby offering a way to reduce considerably the dimensionality of the matrices in high-symmetry space groups. In our treatment here, we ignore internal symmetry.

## 3. Results and discussion

We examined the null modes within the Brillouin zones of all of the known zeolite framework types in their highest symmetry configuration. The goal was to explore the flexibility window of these frameworks and to identify those framework types that have extensive null modes. Not all framework types can be represented as being built from ideal tetrahedra. In several instances, the **MTN** framework being a good example, it is not possible to represent a framework with ideal tetrahedra at maximum symmetry. Thus, the exercise needs to be conducted in the lowest, *P*1, symmetry. In reduced, *P*1, symmetry *all* of the known frameworks exhibit flexibility. Six frameworks, **CZP**, **GOO**, **ISV**, **ITR**, **IWS** and **STW**, however, require two differently sized tetrahedra in order to be flexible. This indicates that the silicate composition of the actual materials inherently needs to be mixed with either aluminium or germanium, which are about 8% larger than the silicon tetrahedra. To accommodate the possibility of mixed compositions, we used the atomic coordinates and unit-cell dimensions provided by the Structure Commission of the International Zeolite Association (IZA) [7]. These coordinates are based on a distance-least-squares minimization against target T–O distances, and T–O–T angles in the maximum symmetry, that is compatible with tetrahedra (some symmetries impose planar tetrahedra). The IZA structures contain nearly ideal tetrahedra, but not perfectly so. We constructed the compatibility matrix to be consistent with the actual *d*_{OO} distances found for each O–O ‘bar’, and we assumed that these imperfect tetrahedra were rigid. We found that 25 of the 206 approved framework types exhibited extensive flexibility, the number of unique null modes growing linearly with the volume of the supercell chosen [27]. The majority of zeolite frameworks exhibit a volume^{2/3} growth rate in number of modes, growing more as a surface term.

In [27], we explored the configurational entropy for the extensive sodalite framework, **SOD**, based on a random sampling of mechanisms for a supercell of the structure, and following them while updating the entity **C**^{−1}**C**^{T}**C****T** as it evolved through the flexibility window. We found that the configurational entropy, computed using equation (2.17), was sharply peaked at the minimum density end of the flexibility window. We argued that the configurational entropy is the highest at the maximum symmetry, which is almost always at maximum unit cell volume (i.e. minimum density). As we show here, the configurational entropy can also vary even if we follow a mechanism while maintaining high symmetry.

The sodalite (**SOD**) framework can be folded while maintaining body-centred cubic symmetry, . The highest symmetry, body-centred cubic , is expressed only at the minimum density, and is violated when any mechanism is followed. Figure 3*a* depicts the **SOD** framework in its lowest density state, with symmetry , and figure 3*b* shows it at its maximum density under symmetry.

The configurational entropy term, *S*(*ρ*), for the **SOD** framework is plotted in figure 4 as a function of density, *ρ*, within the flexibility window (right-hand axis). The symmetry is held as cubic throughout the folding stage of the mechanism. We did this by allowing one degree of freedom in the unit cell, the T–O–T angle. The Si atom is always located at (1/2, 1/4, 0) in cell-relative coordinates, and the unique oxygen site is always at (*x*,*x*,*z*). *x* and *z* are uniquely defined by enforcing O–Si–O to be the tetrahedral angle, and setting Si–O–Si to its desired value, which in turn determines the density. The minimum density is also the maximum symmetry state, cubic , which has the maximum number of degenerate constraints. There is no mechanism that preserves this high symmetry, which occurs only at the lowest density point of the ideal framework. Consequently, the low-density state is also the configuration for the **SOD** framework that has the maximum number of null modes, which is what the configurational entropy term is really measuring. Increasing the density of **SOD** from its minimum density immediately lowers the symmetry, and thus lowers the number of null modes, hence the rapid drop in the configurational entropy.

Also plotted on figure 4 is the T–O–T, intertetrahedral, angle (left-hand axis). There is only one angle possible throughout the structure when the symmetry is maintained as . The T–O–T angle ranges between 160.53° at the lowest density and 114.97° at the highest density. The configurational entropy term is sharply peaked at the minimum density of *ρ*=16.61 T-atoms nm^{−3}. X-ray diffraction studies of pure silica sodalites that contain ethylene glycol show that the framework adopts a density in the range of 17.0–17.4 T-atoms nm^{−3}, with Si–O–Si angles of 159.7–155.2°, in agreement with values we predict for the low-density end of the ideal framework [28,33]. We were unable to locate a structure refinement of pure-silica sodalite free of extra-framework inclusions. This is also the configuration that gives maximum configurational entropy (that is, the maximum number of null modes). However, sodalite is a highly adaptable clathrate framework and its unit cell can vary significantly with composition.

The second example is provided by the **LTA** framework, which also has an extensive count of null modes. Three conformations of a maximum symmetry mechanism () for a single unit cell are depicted in figure 5. Figure 5*a* shows the maximum volume state, at density 13.95 T-atoms nm^{−3}. Upon compression, while maintaining symmetry, the stress is accommodated by rotations of tetrahedra about the oxygen joints. At density 15.88 T-atoms nm^{−3} (figure 5*b*), a significant number of T–O–T angles pass through 180° (arrowed in figure 5*b*). At this density, the framework is in a non-generic configuration, having developed an additional set of degenerate constraints without increasing the space group symmetry of the framework. Upon further densification, the fold continues and the T–O–T angles decrease again (figure 5*c*).

There are three distinct T–O–T angles possible in the **LTA** framework in symmetry . Their values are plotted in figure 6 (left-hand axis). As the mechanism proceeds from the lowest density to higher density, two of the T–O–T angles decrease, and one increases until it reaches 180° at density 15.88 T-atoms nm^{−3}. The configurational entropy (right-hand axis) exhibits a peak at this density. The 180° angles have ‘softened’ (lowered the vibrational frequency) in some of the modes, which raises the configurational entropy at this density. Real pure-silica materials with the **LTA** framework adopt a density of 14.36 T-atoms nm^{−3}, with mean T–O–T angles of 152.4° [34]. Any thermodynamic benefit gained by increasing the configurational entropy at density 15.88 T-atoms nm^{−3} appears to be outweighed by the bonding terms that lower the framework energy. In practice, there is a small framework energy penalty, which our treatment ignores, for frameworks that contain large T–O–T angles. It appears, therefore, that the increase in framework energy induced by T–O–T angles of 180° outweighs any benefit arising from the increase in configurational entropy favoured by that angle.

The third example is offered by the **ABW** framework, which is not an extensive framework. Materials with the **ABW** framework have not been synthesized as a pure silicate, and usually contain a significant concentration of Al atoms. AlO_{4} tetrahedra are about 7.8% larger than SiO_{4} tetrahedra. In **ABW**, the number of null modes grows as volume^{2/3}. When following a mechanism that maintains orthorhombic symmetry, *Imma*, most of the strain is accommodated by a decrease in the length of the unit cell *a* axis (figure 7*a*–*c*). Again, there are three distinct T–O–T angles under this symmetry, one of which is held fixed at 180°. This angle is marked by the asterisks (*) in figure 7*a*,*b*. Similarly to the **LTA** model framework, at intermediate density of 19.38 T-atoms nm^{−3}, the second angle passes close to 180° (marked by the arrows in figure 7*b*). From the **LTA** example, we would expect that the **ABW** model framework is also passing through a non-generic configuration. However, the configurational entropy does not show a peak at this density (figure 8). Inspection of the framework reveals that there are two, non-collinear, folding axes present at the arrowed joints (oxygen atoms). Consequently, the 180° angles that appear at 19.38 T-atoms nm^{−3} do not represent degenerate constraints, and no soft modes evolve. T–O–T angles of 180° do not necessarily generate soft modes.

When the **ABW** framework is started from the lowest density state, which has symmetry *Imma*, but is allowed to follow mechanisms that are constrained to a lower symmetry, *Pna*2_{1}, we find that the structure tends to adhere initially to a higher symmetry *Imma* configuration as the structure densifies. At some point (indicated by the vertical line in figure 8), a bifurcation of modes occurs. The configurational entropy drops rapidly as the density increases. There are now four distinct T–O–T angles, and the window available to the *Pna*2_{1} mechanism is larger, reaching a density of 25.06 T-atoms nm^{−3}. The bifurcation points are not unique, and their appearance possibly depends on symmetry-breaking rounding errors as modes are followed. Other mode-following experiments have reached densities as high as 34.91 T-atoms nm^{−3} (figure 7*d*).

Real materials with the **ABW** framework are observed to adopt a lower symmetry *Pna*2_{1}, and contain no 180° angles. The density adopted is near 19.0 T-atoms nm^{−3}, which is significantly higher than the minimum, 17.41 T-atoms nm^{−3}. As the **ABW** framework does not exhibit an extensivity in null mode count, there is no configurational entropic advantage to adopting the lowest density. As the plot in figure 8 shows, the configurational entropy is essentially flat over a wide density range. Further, the framework is sufficiently underconstrained that the symmetry can remain maximal over this range. It appears that the framework energy reduction, obtained by avoiding high-energy 180° Si–O–Si angles, outweighs any entropic advantage of maintaining the higher symmetry. However, it is also observed that many silicate structures tend to go to higher symmetry, or higher configurational entropy, at high temperature. A well-known example is the low quartz to high quartz transition [35], where thermal vibrations probably conform to rigid-body tetrahedral rotations [36] and explore higher entropy (i.e. higher symmetry) states.

## 4. Conclusion

By treating zeolite frameworks as a mechanical periodic truss, with rigid bars connecting the oxygen joints, we have been able to apply rigidity theory to explore the density of mechanisms within the flexibility window of the structure. We find that about 12% of the known frameworks exhibit an extensive growth of mechanisms, where the number of independent mechanisms grows linearly with the volume of the supercell. Such frameworks potentially have a thermodynamic advantage at elevated temperatures in that the ‘floppy’ mechanisms may contribute a significant term to the configurational entropy of the framework. We have merged the methods of rigidity theory with the phonon theory of solids at high temperatures to obtain an expression that approximates the configurational entropy of a flexible periodic truss.

All of the known zeolite frameworks exhibit a flexibility window, even though most frameworks do not exhibit an extensive flexibility (i.e. the number of mechanisms grows slower than the volume of the supercell, typically as *V* ^{2/3} or *V* ^{1/3}). By comparison with experimental results, we conclude that the configurational entropy in already formed zeolite crystals is likely to be only a small contribution to the framework's free energy. We conclude that the T–O–T bending forces are more important in determining the equilibrium configuration of zeolite frameworks than is the configurational entropy. The framework free energy is dominated by electrostatic terms, which ensure that the SiO_{4} tetrahedra are close to ideal. However, at elevated temperatures, and possibly during nucleation, the configurational entropy, associated with the floppy modes arising from the flexibility of the tetrahedral network about the T–O–T linkages, may play a more important role in determining the final structure.

All zeolite frameworks support at least a volume^{2/3} growth of null modes as a function of supercell volume. This suggests that during nucleation and growth, where the growing particle's surface is important, such surface-extensive modes could still provide a significant thermodynamic advantage. Zeolite frameworks are metastable to quartz, and therefore do not represent an absolute thermodynamic minimum. Zeolite growth is governed by kinetics. The flexibility of nuclei may provide a thermodynamic advantage at the nucleation stage, where many different shapes and conformations can be adopted by a single topology. It is not the configurational entropy of the final crystal, but that of the nuclei, that may be important in zeolite synthesis. The role of flexibility in the nucleation of small clusters merits closer inspection.

An important conclusion from these studies is that the T–O–T energy terms seem to be more important in determining equilibrium structures of bulk zeolites than we originally appreciated. That is, the forces at the spherical joints need to be included in any analysis of the energetics. However, the thermodynamic advantages of flexibility may be significant for stabilizing nuclei during the early stages of synthesis.

## Funding statements

We are grateful for financial support from the National Science Foundation grant no. DMR-0835605. I.R. thanks the Institute for Advanced Study and the Institut für Mathematik, Techniche Universität Berlin, for support.

## Acknowledgements

Some of the ideas explored in this paper were inspired in part by discussions with Simon Guest at the workshop ‘Geometric constraints with applications in CAD and biology’, organized by Ileana Streanu at Bellairs Research Institute, Barbados, January 2009.

## Appendix A. Handling the degrees of freedom introduced by changes in the unit cell.

Guest & Hutchinson [25] present explicit equations showing the relationship between the 3*N*+6 elements of **C** when the unit cell parameters change. A compact notation can be developed if we use the crystallographic *relative* coordinates, as opposed to the Cartesian coordinates. The relative coordinates *r*_{r}=(*x*_{r},*y*_{r},*z*_{r}) are the fractional coordinates along each of the axes of the generally oblique unit cell. Thus, (1/2, 1/2, 1/2) is a point at the centre of the cell. *r*_{r} can be related to the Cartesian coordinates, *r*_{c}=(*x*,*y*,*z*), used earlier, via a matrix transformation of the type
A1

The choice of transformation is not unique. There is an infinite number of possible transformation matrices ** A**. A useful one is the transformation that aligns the

*z*-axis with

**, and which places the**

*c**y*-axis in the

**–**

*b***plane, thus committing the**

*c**x*-axis to be perpendicular to the

**–**

*b***plane, A2where A3Usually, is found from , although a separate expression exists. The inverse of**

*c***exists provided the cell volume is not zero.**

*A*For small changes in the relative coordinates *δ**r*_{r}, and for small changes in the unit cell parameters *a*, *b*, *c*, *α*, *β*, *γ*, the change in Cartesian coordinates is
A4 is the change in Cartesian coordinates owing to the movement of the atoms relative to the periodic ‘medium’. The second term is the additional ‘inflationary’ shift owing to the change in the unit cell metric. This latter equation can be expressed more explicitly as
A5All reference to the relative coordinates *r*_{r} has been removed. It should not be overlooked that the periodic displacements *X*_{ij}, etc., in equation (2.2) also have non-zero components *δX*_{ij}, etc., and should be included. The relations for these are found from equation (A 5) by setting and *r*_{c}=(1,0,0), or (−1,0,0) as appropriate for each bar, for *X*_{ij}, etc.

## Footnotes

One contribution of 14 to a Theo Murphy Meeting Issue ‘Rigidity of periodic and symmetric structures in nature and engineering’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.