# A Biased Monte Carlo Scheme for Zeolite Structure Solution

## I Introduction

Zeolites continue to be synthesized at a furious pace. Crucial to the development of the field of zeolite science is the ability to determine the structure of newly-synthesized materials: Structure is sought after not only to understand the performance of newly synthesized catalysts but also to propose rational syntheses of homologous materials with tailored performance. Roughly 118 framework structures have been reported, yet another several dozen distinct synthetic zeolites remain unsolved in the patent literature. Perhaps soon the techniques of diversity synthesis will be introduced to the field, with a tremendous explosion in the number of new, unsolved synthetic zeolites.

Zeolites are crystalline microporous materials that have found a wide range of uses in industrial applications. They are used as catalysts, molecular sieves, and ion-exchangers. A typical example is ZSM-5, shown in Fig. 1, and used as a cracking co-catalyst in the refinement of crude oil. The pore structure of this particular zeolite is 3-dimensional, as is that of the zeolites chabazite, A, X, Y, and beta, but there are zeolites with 1-dimensional pores, such as cancrinite, zeolite L, and AlPO<sub>4</sub>-5, and 2-dimensional pores, such as deca-dodecasil 3R, TMA-E(AB), and heulandite. Classical zeolites are aluminosilicates. The basic building block is a TO<sub>4</sub> tetrahedron. Usually T = Si, although substitution of the silicon with aluminum, phosphorus, or other metals is common. The tetrahedral species is commonly denoted by T when one is concerned with structural, rather than chemical, properties of the zeolite.

The derivation of an atomic-scale model of the framework crystal structure of a newly-synthesized zeolite is a non-trivial task. The difficulty stems primarily from the polycrystalline nature of most zeolite samples, with crystallite sizes typically below 5 $`\mu `$m. Notable improvements in single-crystal diffraction techniques have been made, primarily through the use of synchrotron X-rays, but fundamental limitations still exist. A typical powder pattern is shown in Fig. 2. Most zeolites have been solved to date through physical model building efforts. A limited number of zeolites have been solved with conventional crystallographic methods (e.g., see Ref.  for a complex example). These methods attempt to deconvolute the powder diffraction data into unique reflections and then apply direct methods to determine a structure. This approach requires a high quality scattering dataset for success, a dataset that usually is not obtainable from zeolite samples. Very few zeolite structures have been solved via conventional crystallographic methods, which indicates the limited applicability of this approach. Recently, a few zeolites have been solved with electron diffraction methods, but this is a laborious approach.

A direct, real-space method for zeolite structure solution from powder diffraction data has been proposed. The advantage of this method is that it requires only easily available powder data, and it incorporates little preconceived bias of the investigator in the structure solution scheme. In this approach, simulated annealing with a simple Metropolis perturbation step is used to minimize a zeolite figure of merit. This approach has a nearly 90% success rate on materials with 6 or fewer crystallographically-distinct T-atoms. To date, at least four other groups have used this approach to solve new zeolite structures. On more complex structures, however, the method becomes unwieldy, generating many hypothetical structures and often failing to find the correct one.

In this paper, we apply powerful new ideas from Monte Carlo to the problem of zeolite structure solution via powder diffraction data. The key step is to define a cost function that is a function of the atomic positions within the crystalline unit cell and that is minimized by the structure corresponding to the experimental material. Structure solution from powder data is a fundamentally challenging problem, due to the presence of many local minima and large barriers in the cost function. We show that a combination of simulated annealing and biased Monte Carlo is able to overcome the barriers in most cases. In the most difficult cases, we show that the new method of parallel tempering is superior to simulated annealing and is able to overcome the barriers. This paper is organized as follows. In Section II we introduce the figure of merit for an unknown zeolite sample. In Section III we discuss the Monte Carlo sampling of this figure of merit. In Section IV we show the results of applying our method to all publicly known zeolites. In Section V we discuss the results of the method and possible extensions. We draw our conclusions in Section VI.

## II The Figure of Merit

For a given zeolite sample with known unit cell size, cell parameters, symmetry, and density we want to construct a figure of merit. By definition, the global minimum of the figure of merit should correspond to the structure of the zeolite sample that we are investigating. The $`n_{\mathrm{unique}}`$ unique atoms are placed in the unit cell. Each of the $`n_{\mathrm{symm}}`$ symmetry operators generates an atom position when applied to a unique T-atom, and so there are at most $`n_{\mathrm{unique}}\times n_{\mathrm{symm}}`$ atom positions in the unit cell. The figure of merit is a function of all of these parameters. We keep fixed the cell parameters, the space group symmetry, and the number of unique T-atoms: the only variables are the positions of the T-atoms in the unit cell. The figure of merit is defined as

$`H`$ $`=`$ $`\alpha _{\mathrm{T}\mathrm{T}}H_{\mathrm{T}\mathrm{T}}+\alpha _{\mathrm{T}\mathrm{T}\mathrm{T}}H_{\mathrm{T}\mathrm{T}\mathrm{T}}`$ (1)   
$`+`$ $`\alpha _{\mathrm{T}\mathrm{T}\mathrm{T}}H_{\mathrm{T}\mathrm{T}\mathrm{T}}+\alpha _\mathrm{D}H_\mathrm{D}+\alpha _\mathrm{M}H_\mathrm{M}+\alpha _\mathrm{L}H_\mathrm{L}`$ (2)   
$`+`$ $`\alpha _{\mathrm{NB}}H_{\mathrm{NB}}+\alpha _{\mathrm{PXD}}H_{\mathrm{PXD}}+\alpha _{\mathrm{PND}}H_{\mathrm{PND}}.`$ (3)   
The figure of merit is defined for any arrangement of T-atoms, even ones that are far from resembling a zeolite. The lower the value of $`H`$, the more the structure resembles a zeolite of the given cell size, symmetry, and density. In general, one starts from a random configuration of T-atoms and seeks the minimum of $`H`$ by moving the T-atoms suitably. Each term in $`H`$ represents a particular contribution, and the $`\alpha _i`$ are the relative weights. These weights are optimized according to the rate of success of the method on one trial structure and are generally kept fixed in the following.

It must be stressed that $`H`$ is not the thermodynamic energy. For one, we neglect the bridging oxygens and any cations or adsorbed molecules that may be present in the real zeolite. Furthermore, we will see that the geometric terms of $`H`$ do not describe interactions or pseudo-interactions, but are constructed to reproduce a distribution of distances and angles observed in known zeolites. We enforce the space group symmetry, and so symmetry-related T-atoms move collectively. In addition, we enforce crystalline order so that we can limit our description of the material to a single unit cell.

One may wonder whether the explicit inclusion of bridging oxygens, or whether a more detailed, quantum mechanical description of the system, may be needed for an accurate description of the zeolite. In fact, Eq. (3) captures the relevant features necessary to describe and to predict the framework structure of the zeolite crystal. Once the positions of the T-atoms are available, one can use more refined methods and more detailed models to refine the structure or to study other properties that may be of interest.

The different contributions to Eq. (3) are of three types: the geometric terms, the density terms, and the diffraction terms. The first three terms of $`H`$ are geometric, and they are obtained by histogramming the T–T distances and the T–T–T angles of 32 known high silica zeolites. The T–T distance is sharply distributed around 3.1 Å, and the T–T–T angles are distributed around 109.5<sup>o</sup>, as one would expect for a tetrahedrally coordinated species. The angles and distances are shown in Fig. 3. The angle $``$T–T–T$``$ is the average of all of the angles around a given T-atom. We are interested in sampling configuration space according to Boltzmann statistics. This leads us to define potential energies that reproduce the observed histograms when sampled according to $`\mathrm{exp}(\beta H)`$ at a particular value of the inverse temperature $`\beta =1/T`$. For simplicity, we set Boltzmann’s constant to unity. In practice we take the logarithm of the histogrammed data and fit a spline through them. The potential energies are then extended to ranges of angles and distances beyond the observed ones. The resulting curves are shown in Fig. 4a, 4b, and 4c. The correct histograms are obtained only when the distribution is sampled at the inverse temperature $`\beta `$ used for the inversion. This value is arbitrary and simply defines what low temperature means.

The $`H_{\mathrm{NB}}`$ and $`H_\mathrm{L}`$ terms account for the 4-connectedness of silicates. These terms are defined to be non-zero and positive whenever a T-atom happens to have greater than or fewer than 4 first neighbors respectively, where a neighbor is defined as an atom closer than 4 Å. If the number of neighbors is fewer than 4, we simply assign a progressively larger weight to the atom. Following Ref.  we use the values in Table I. The case in which there are more than 4 neighbors is treated differently. We define a list of the $`N`$ neighboring atoms. We choose 4 of them and treat them as connected neighbors, adding the contributions from the potential energies discussed above. We include a repulsive potential energy, Fig. 4d, for the remaining unconnected neighbors. Since it is important to choose the best connected neighbors for a configuration, we exhaustively search for the combination of 4 connected and $`N4`$ unconnected atoms that minimizes the energy associated with the central $`T`$ atom. In our experience, this procedure gives better results than does simply picking the neighbors according to their distance or bonding all neighbors indiscriminately and assigning an extra weight to over-coordinated atoms.

The $`H_\mathrm{M}`$ term favors merging. Merging occurs in crystals whenever a particular atom sits on a special position, a position invariant under one or more symmetry operations other than the identity. Since our method assigns positions in a stochastic way, it is unlikely that we would find an atom exactly on a special position. Therefore, we define a merging range with a typical value of $`r_\mathrm{M}=0.8`$ Å. Two or more symmetry-related atoms that fall within this distance are merged. When this condition occurs, all the atoms that fall within the range are replaced by a single atom at the position of their center of mass. The merged position is not a new position for the system; it is merely the position used for the calculation of the figure of merit. $`H_\mathrm{M}`$ gives a negative, favorable energy to merged atoms. This energy is linearly proportional to the merging distance, i.e. the distance between the original and the merged position. Merging must be favored since the figure of merit is proportional to the number of atoms in the unit cell, and placing an atom on a special position lowers the number of actual atoms in the unit cell and will typically remove a favorable energy contribution. Merging is necessary whenever the number of T-atoms derived from the experimental density, $`n_\mathrm{o}`$, is fewer than the number created by symmetry, $`n_{\mathrm{unique}}\times n_{\mathrm{symm}}`$. Merging is disallowed whenever $`n_\mathrm{o}=n_{\mathrm{unique}}\times n_{\mathrm{symm}}`$. If merging is allowed, the number of T-atoms in the unit cell can change. In fact since the contribution $`H_\mathrm{M}`$ that favors merging is negative, there are situations in which some atoms collapse on a highly symmetric special position, lowering the density far below the observed one. In order to enforce the observed density we include the term

$$H_\mathrm{D}=(n_\mathrm{T}n_\mathrm{o})^2,$$
(4)   
where $`n_\mathrm{T}`$ is the actual number of atoms in the unit cell after merging.

We now discuss some technical details of the figure of merit thus defined. Symmetry can be exploited in the computation of the figure of merit in several ways. Symmetry-related atoms will have the same geometric structure surrounding them, so the energy is simply $`n_{\mathrm{symm}}`$ times the energy of a single atom, unless merging has occurred. In any case, the energy is proportional to the actual number of atoms in the unit cell. To speed up the calculation of the energy, we divide the unit cell with a grid. Each atom belongs to a box in the grid, and the grid is designed so that first neighbors to an atom will be in the same box as the atom or in one of the 26 neighboring boxes. The computational effort to calculate the energy is reduced from $`O(n_{\mathrm{unique}}^2)`$ to $`O(n_{\mathrm{unique}})`$ through use of the grid. The grid must be updated every time an atom is moved, but since only one unique atom is moved at one time, we create $`n_{\mathrm{unique}}`$ distinct grids. This allows us to update only the grid associated with a moved atom, leaving the others unchanged and further reducing the computational effort used per attempted move. The figure of merit we defined is highly non-local. Symmetry operations applied to the position of the unique atom can, in principle, generate atoms anywhere in the unit cell. This means that the energy for all of the unique atoms must be recalculated every time a single unique atom is moved.

The diffraction terms, $`H_{\mathrm{PXD}}`$ and $`H_{\mathrm{PND}}`$, incorporate in the figure of merit experimental information that may be available about the zeolite. A typical PXD pattern is shown in Fig. 2. Let us assume that a skilled crystallographer has collected high resolution powder X-ray scattering data on a zeolite powder sample and has succeeded in indexing the resulting pattern. We then have available a list of Bragg reflections with Miller indices $`(hkl)`$ and relative intensities. For a given arrangement of atoms in our model unit cell, we can compute the relative intensities for the same list of reflections using standard formulas. The intensity of a reflection, in arbitrary units, is given by

$$I(hkl)=p(\theta )|F_{hkl}|^2,$$
(5)   
where $`2\theta `$ is the angle of the Bragg reflection, and the first factor is the polarization term. The polarization term is $`p(\theta )=[1+\mathrm{cos}^2(2\theta )]/[2\mathrm{sin}(\theta )\mathrm{sin}(2\theta )]`$ for X-rays and $`p(\theta )=1/[2\mathrm{sin}(\theta )\mathrm{sin}(2\theta )]`$ for neutrons. The scattering amplitude is

$$F_{hkl}=\underset{j=1}{\overset{n_\mathrm{T}}{}}f_j(𝐤)o_j\mathrm{exp}(B_j𝐤^2/4)\mathrm{exp}(i𝐤𝐱_j),$$
(6)   
where

$`𝐤`$ $`=`$ $`h𝐛_1+k𝐛_2+l𝐛_3`$ (7)   
$`𝐱_j`$ $`=`$ $`m_j^{(1)}𝐚_1+m_j^{(2)}𝐚_2+m_j^{(3)}𝐚_3.`$ (8)   
Here, the $`𝐚_i`$ are the crystal axes, the $`m_j^{(i)}`$ are the crystallographic coordinates, and the $`𝐛_i`$ are the reciprocal lattice vectors. The $`f_j(𝐤)`$ are the form factors for the given atomic species; the $`o_j`$ are the occupancy numbers, which account for cell positions not always filled with an atom or filled with atoms of different type with different probabilities; and the $`B_j`$ are the isotropic Debye-Waller factors that account for thermal vibrations in the lattice.

Our description of the contents of the cell is approximate, since the oxygens and cations are excluded. Also, since the structure of the crystal is unknown, we do not have information about the occupancies or Debye-Waller factors. Therefore, we set $`o_j=1`$ and $`B_j=1/2`$ in all our trials. These limitations imply that our calculated diffraction pattern cannot exactly match the observed one even if we locate perfectly all of the framework T-atoms. Nonetheless, we can capture the relevant features of the pattern. The contribution of the oxygens to the diffraction pattern is less important than that of the silicons, since $`f_j`$ is roughly proportional to the atomic number. Indeed, it proves better to leave the oxygens out than to include them at the midpoints between the T-atoms. The contribution of non-framework species to the reflection intensities is suppressed since they usually have large Debye-Waller factors, being more loosely bound than the framework species, and often have fractional occupancies.

The presence of multiple reflections at angles closer than the resolution obtainable even with the best synchrotron radiation sources is one of the major challenges to the use of powder data. This occurrence is common for zeolite samples. In order to compare the computed intensities with the experimental ones, we define a composite peak at the average angle and place all of the intensity of the multiple reflections into this composite peak.

To measure how well a particular configuration of T-atom positions can match the experimental powder pattern, we define the quantity

$$H_{\mathrm{PXD}}=\frac{1}{N}\underset{s}{\mathrm{min}}\left[\frac{_i(I_i^{\mathrm{obs}}sI_i^{\mathrm{calc}})^2/\omega _i}{_i1/\omega _i}\right],$$
(9)   
where $`i`$ runs over all of the $`N`$ peaks, which may be composites, the $`\omega _i`$ are the weights, and $`s`$ is a global scaling factor. We make a similar definition for any available neutron data. The intensities are relative, and we scale the experimental intensities so that the largest one is 1000. The explicit expression for the global scaling factor, $`s_{\mathrm{min}}`$, is

$$s_{\mathrm{min}}=\frac{_i(I_i^{\mathrm{obs}}I_i^{\mathrm{calc}})/\omega _i}{_i\left(I_i^{\mathrm{calc}}\right)^2/\omega _i}.$$
(10)   
The weights $`\omega _i`$ are associated to each peak according to the following criterion: if the scaled intensity is less than 90, $`\omega _i=1`$; if $`90<I_i150`$, $`\omega _i=2`$; if $`150<I_i300`$, $`\omega _i=3`$; and $`\omega _i=4`$ otherwise. These weights account for uncertainty both in the data and our fit to the data, both roughly proportional to the intensity itself.

The numerical implementation of the calculation of the diffraction terms can exploit the presence of the space group symmetry. Since we move one unique T-atom at a time, only the contribution of the unique T-atom and of all of its symmetric images will change in $`H_{\mathrm{PXD}}`$ and $`H_{\mathrm{PND}}`$. To calculate energy changes, only these terms need to be reevaluated. We can express the sum Eq. (6) in the following way:

$$F_{hkl}=\underset{\gamma }{}F_{hkl}^{(\gamma )},F_{hkl}^{(\gamma )}=\underset{j\gamma }{}f_j(𝐤)o_j\mathrm{exp}(B_j𝐤^2/4)\mathrm{exp}(i𝐤𝐱_j).$$
(11)   
Only one of the terms of the first sum changes when a single unique T-atom is moved.

A good choice of the weights $`\alpha _i`$ in Eq. (3) is crucial for the success of method. The diffraction terms are very sensitive to the positions of the atoms. In other words, they make for a very rough energy profile. Even a small change in the position of one atom and its symmetric images can change $`H_{\mathrm{PXD}}`$ and $`H_{\mathrm{PND}}`$ substantially. The geometric potentials of Fig. 4 are quite smooth. The density terms are also relatively smooth, even when the move involves a change in the merging. We want to mix the geometric terms with the diffraction terms in such a way that the roughness of the diffraction terms is smoothed out. In order to find the best values of these weights, we tried several different combinations on one trial structure. It was reassuring to find that the rate of success was not sensitive to small changes in the parameters. A good choice for the $`\alpha _{\mathrm{PXD}}`$ or $`\alpha _{\mathrm{PND}}`$ weight is between 1 and 2. The diffraction terms are extensive, in the sense that the number of reflections is roughly proportional to the total number of atoms, as are the other terms in the figure of merit, but in few cases we found it necessary to increase or decrease their importance with respect to the other terms. The values of these weights are $`\alpha _{\mathrm{T}\mathrm{T}}=1`$, $`\alpha _{\mathrm{T}\mathrm{T}\mathrm{T}}=1`$, $`\alpha _{\mathrm{T}\mathrm{T}\mathrm{T}}=2.0`$, $`\alpha _{\mathrm{NB}}=1.5`$, and $`\alpha _\mathrm{L}=1`$. The merging term $`H_\mathrm{M}`$ was in almost all cases fixed to be zero at $`r_\mathrm{M}=0.8`$ Å and $`300`$ at $`r_\mathrm{M}=0`$ Å, and the associated weight was set at $`\alpha _\mathrm{M}=1`$. The density weight $`\alpha _\mathrm{D}`$ was usually around $`30`$, in order to ensure that the proper density was reached at low temperatures.

How the figure of merit depends on the atom coordinates is shown in Fig. 5. The curve corresponds to a one dimensional slice of the full profile, obtained by sliding one coordinate of one T-atom across the unit cell of the zeolite faujasite. The profile is very rough, with many narrow valleys. It is crucial to notice that the position of the atom in the faujasite structure corresponds to the global minimum of the curve, at roughly $`m_1^{(1)}0.6`$. If we were solving this material, this position would be the minimum that we would have to locate in order to solve the structure. Clearly the figure of merit is very rough, and we will need a powerful simulation protocol to perform the many-dimensional global optimization.

The figure of merit defined by Eq. (3) may possess invariance under translations of the unit cell in particular directions. This reflects the fact that for some space group symmetries the position of the unit call may be arbitrary. This is clear, for instance, in the case of the space group P1. In this space group, the unit cell can be moved freely in all directions, without changing the description of the system. Our figure of merit would be invariant under a simultaneous, arbitrary, and continuous translation of the atoms. The space groups that have this type of freedom are called polar, and the directions in which the unit cell can be moved are the polar directions. In other instances, the atoms within the unit cell can be moved a discrete amount in one direction, such one-half the unit cell, and still lead to the same crystal structure. In both of these cases it is useful to eliminate these translationally invariant modes. It is straightforward to identify a polar group and the polar directions and to define a projection operator that will restrict proposed moves to an orthogonal subspace of the polar directions. It is enough to select one atom that can move only along the orthogonal subspace to break the polar symmetry. In the case of unit cells with two or more choices of cell setting, such as occurs in the framework LTL, the ambiguity can be eliminated by choosing an appropriate asymmetric unit and limiting the movement of the unique T-atom to that cell. Use of an asymmetric unit, however, was found not to be necessary for the success of the method.

## III The Monte Carlo

In this section we will describe in some detail the Monte Carlo algorithm that we use to sample the figure of merit. We will first make some general comments about biased Monte Carlo importance sampling. We then describe simulated annealing and parallel tempering.

Monte Carlo methods have been used extensively since their inception to sample equilibrium probability distributions of systems with many degrees of freedom. The key step is to define a Markov process that evolves the system from configuration to configuration. As long as this Markov process satisfies certain properties, one is assured that after $`N`$ steps, time averages will approximate ensemble averages to within a relative error of $`1/\sqrt{N}`$. Specifically, if the Markov process is ergodic and regular and satisfies detailed balance, it can be shown that the limiting probability distribution is the one we seek. The proof uses the Perron-Frobenius theorem and the fact that a matrix obeying detailed balance has a complete set of eigenvectors (see, for example, Ref. ). A more general proof shows that the method need satisfy only the weaker balance condition.

One of the shortcomings of the traditional Metropolis method is that it does not use any information about the energy landscape around the current configuration when picking trial moves. Oftentimes, the proposed move brings the system to regions of configuration space that are high in energy, and the move is rejected. These rejected moves hinder effective sampling of the Boltzmann distribution.

Biased Monte Carlo methods have been shown to improve sampling in many cases. They were originally introduced to lead to more efficient simulations of complex liquids. The basic idea is to probe the configurations around the current one and to propose moves that are more likely to be accepted. In our case the biased move proceeds as follows. Let us call the current configuration $`A_1`$. We extract $`k`$ random displacements, $`\mathrm{\Delta }𝐱_i`$, which define $`k`$ proposed new configurations, $`B_i`$. These moves are extracted from a Gaussian distribution

$$p_i^{\mathrm{int}}=\frac{\mathrm{exp}[\mathrm{\Delta }𝐱_i^2/(2\sigma ^2)]}{[2\pi \sigma ^2]^{3/2}}.$$
(12)   
We construct the Rosenbluth weight $`W`$, defined as

$$W(n)=\underset{i=1}{\overset{k}{}}\mathrm{exp}[\beta H(B_i)],$$
(13)   
and we assign a normalized probability

$$p_i^{\mathrm{ext}}=\mathrm{exp}[\beta H(B_i)]/W(n)$$
(14)   
to each configuration $`B_i`$ (see Fig. 6). We randomly select one of these configurations, $`B_n`$, according to its probability. The configuration $`B_n`$ is our proposed move. Clearly the lower the energy is, the more likely the configuration will be selected. In order to satisfy detailed balance, we must modify the acceptance probability of the proposed move. This requires us to calculate the likelihood of the reverse move $`B_nA_1`$. The super detailed balance condition, which ensures detailed balance, can be satisfied by defining a set of $`k1`$ new trial moves, $`A_j`$, from the proposed configuration, $`B_n`$. The set $`\{A_1,A_j\}`$ defines the reverse Rosenbluth weight

$$W(o)=\mathrm{exp}[\beta H(A_1)]+\underset{j=2}{\overset{k}{}}\mathrm{exp}[\beta H(A_j)],$$
(15)   
and the normalized probability of selecting the reverse move is

$$p_o^{\mathrm{ext}}=\mathrm{exp}[\beta H(A_1)]/W(o).$$
(16)   
The super detailed balance condition can now be written as

$$\pi (A_1)T(A_1B_n)\mathrm{acc}(A_1B_n)=\pi (B_n)T(B_nA_1)\mathrm{acc}(B_nA_1),$$
(17)   
where $`\pi (A)\mathrm{exp}[\beta H(A)]`$ is the limiting distribution that we want to sample, and the probability of accepting the proposed move is acc$`(A_1B_n)`$. The forward transition probability $`T(AB_n)`$ is just the probability of selecting the configuration, $`p_n^{\mathrm{int}}p_n^{\mathrm{ext}}`$, and the reverse transition probability is $`p_o^{\mathrm{int}}p_o^{\mathrm{ext}}`$. The super detailed balance condition is, then,

$$\frac{\mathrm{acc}(A_1B_n)}{\mathrm{acc}(B_nA_1)}=\frac{\pi (B_n)}{\pi (A_1)}\frac{p_o^{\mathrm{int}}}{p_n^{\mathrm{int}}}\frac{p_o^{\mathrm{ext}}}{p_n^{\mathrm{ext}}}=\frac{W(n)}{W(o)}.$$
(18)   
A reasonable choice for the acceptance probability is

$$\mathrm{acc}(A_1B_n)=\mathrm{min}(1,\frac{W(n)}{W(o)}).$$
(19)   
This class of biased moves significantly improves the sampling of our algorithm with respect to the simple Metropolis scheme. We found that a biased move with $`k=5`$ works well.

The figure of merit described in the previous section gives a quantitative measure of how well a particular arrangement of atoms resembles a zeolite. We are interested in minimizing the figure of merit in order to find the most reasonable arrangements. Biased Monte Carlo alone is unable to sample efficiently the rough figure of merit at low temperatures. Sampling can often be achieved, however, with simulated annealing. In this approach, a series of simulations at progressively lower temperatures is performed, and the distribution at each temperature is sampled using a Monte Carlo method, with or without biasing. The simulation is started from a high temperature and the temperature is progressively reduced according to an annealing scheme. In general the temperature is kept unchanged for a fixed number, $`N`$, of Monte Carlo steps. After these $`N`$ steps the temperature is reduced according to $`T^{}=\kappa T`$, where $`\kappa <1`$, and this cycle is repeated until the temperature is such that most Monte Carlo moves are rejected, and the system is effectively frozen.

The width of the distribution of proposed moves, $`\sigma `$ in Eq. (12), can be adjusted during the simulated annealing run. In general, for a fixed value of the temperature one chooses $`\sigma `$ so that a reasonable number of moves are accepted. We call the ratio of accepted to attempted moves the acceptance ratio $`g`$. On the one hand, if the trial moves are small, then most moves are accepted, $`g1`$, but the effect on the energy is minimal, and the probability distribution is sampled ineffectively. On the other hand, if the trial moves are large, then most moves are rejected, and the sampling is also ineffective. In simulated annealing the issue of acceptance is further complicated by the fact that the temperature is changed. We found it convenient to fix a target acceptance ratio, $`g_\mathrm{t}`$, and to adjust the size of the proposed move distribution, $`\sigma `$, so as to make the actual $`gg_\mathrm{t}`$. In general we start with a large width, $`\sigma =3`$ Å, and we lower it during the annealing to values around $`\sigma 0.51.0`$ Å at low temperature. We use a proportional control scheme

$$\sigma ^{}=\sigma [1+ϵ(gg_\mathrm{t})],$$
(20)   
to adjust $`\sigma `$ each time we reduce the temperature. Even though $`\sigma `$ lags with respect to the temperature, since we use the $`g`$ measured at the higher temperature, we found this control scheme to be effective.

The initial high temperature is found by fixing $`\sigma `$ and performing short trial runs with a Metropolis Monte Carlo. We always start from a completely random initial condition. The temperature is doubled until the fraction of accepted moves during the trial run exceeds a given threshold value. We found that a threshold of $`g=0.5`$ is always sufficient to get to a high enough temperature. Once this initial stage is completed, we thermalize the structure at this fixed temperature with the biased moves. This ensures that we lose track of the initial condition. We then start cooling the system according to a preset annealing schedule. In all but a few cases we used $`N=200`$ and $`\kappa =0.8`$. A typical annealing energy trace is shown in Fig. 7.

We will show in Section IV how effective the combination of simulated annealing and biased moves is in finding the correct frameworks of known zeolites. Typically one or a few runs at most are needed to solve a structure. If the first run is not successful, we try again with different initial positions and random seed, and eventually the correct structure is found. For complex structures with many unique atoms, $`n_{\mathrm{unique}}8`$, this approach sometimes fails to converge to the correct structure within a reasonable time. Of course, one could try to use different annealing schedules, or try with a greater number of different initial conditions. In principle, nothing prevents this method from finding the correct solution. Nonetheless, slow equilibration is more than just a technical detail. For one, during simulated annealing the system is not at equilibrium, since the temperature is reduced at regular intervals. For another, once the system falls in a local minimum in the rough energy profile (Fig. 5), and the temperature is too low for the system to escape via fluctuations in a finite number of steps, the system is stuck.

One other method stands out as a good candidate for sampling probability distributions with complicated landscapes: parallel tempering. This method was developed as an effective Monte Carlo procedure for the study of systems with large free energy barriers. This method was later applied to spin glasses, self-avoiding random walks, lattice QCD, and studies of biological molecules. Following Ref. , we call the method parallel tempering, for its similarity to simulated tempering, a related method also proposed as an improvement on simulated annealing. J-walking is a similar method, often used in molecular optimization problems. J-walking is not an exact Monte Carlo scheme, due to the non-Markovian reuse of configurations. While none of the implementations of parallel tempering cited above is Markovian at the level of a single move, it is a simple matter to make such an implementation. Indeed, taking care with the definitions allows one to understand how to optimize the parallel tempering method by more frequently updating the systems at lower temperature and with longer autocorrelation times.

The idea of parallel tempering is to consider $`n`$ systems, each in a canonical ensemble, and each at a different temperature. We define the instantaneous configuration of system $`i`$ at Monte Carlo step $`t`$ to be $`C_i(t)`$. Each system $`i`$ has a different temperature $`T_1<T_2<\mathrm{}<T_n`$, where $`T_1`$ is the low temperature that we want to sample, and $`T_2,\mathrm{},T_n`$ are higher temperature systems that aid in the sampling. The extended canonical ensemble is given by

$$𝒬=\underset{i=1}{\overset{n}{}}𝒬_i,$$
(21)   
where $`𝒬_i`$ is the canonical partition function,

$$𝒬_i=\underset{\{C_i\}}{}\mathrm{exp}[\beta _iH(C_i)].$$
(22)   
We introduce a swap move, which proposes the exchange of two copies at different temperatures. The proposed move is accepted according to the Metropolis rule. We compute the action difference that the swap move introduces,

$$\mathrm{\Delta }S=\beta _jH(C_i)+\beta _iH(C_j)\beta _jH(C_j)\beta _iH(C_i)=(\beta _j\beta _i)[H(C_i)H(C_j)].$$
(23)   
To satisfy detailed balance, we accept the move with probability

$$p=\mathrm{min}[1,\mathrm{exp}(\mathrm{\Delta }S)].$$
(24)   
Typically we consider swaps between adjacent temperatures, $`j=i+1`$. For a good choice of temperatures, swaps will be accepted with a significant probability. We show in Fig. 8 a schematic drawing of the swapping process.

Parallel tempering allows the system to escape local minima by swapping with the systems at higher temperature. The choice of temperatures should be such that the high temperature, $`T_n`$, is great enough so that the extended ensemble can effectively surmount the free energy barriers. The intermediate temperatures create a ladder that the system uses to climb over the barriers. It is important to notice that the extended ensemble is precisely defined by Eq. (21). This means, for example, that system $`i`$ samples the canonical ensemble at temperature $`T_i`$. We satisfy detailed balance because of Eq. (24). The displacement and swapping moves are clearly ergodic, in principle. If our moves are defined so as to produce a Markov process, then we are guaranteed to sample the extended ensemble in Eq. (21).

When there is more than one kind of update rule in a Monte Carlo simulation, the moves must be selected randomly in order to have a Markov chain on the level of a single move. Of course, one is free to pick the relative probabilities of selecting each type of move. Our implementation of parallel tempering selects the moves at random. We start by selecting one of the systems at random. We then randomly decide whether to make a swap move or a displacement move. We have found that choosing a displacement move 90% of the time leads to efficient sampling. If a displacement move is selected, one of the unique atoms, chosen at random, is updated. If a swap move is selected, we attempt to swap the chosen system with the system at higher temperature. Since the systems at low temperature are slower to evolve under the Monte Carlo sampling, we pick these systems more frequently than the ones at higher temperature. We typically allow the two lowest temperatures to be updated twice as frequently, leading to more swapping and more displacement moves for these systems. In the general case, the update frequencies should be proportional to the autocorrelation times of the respective systems, as measured in the parallel tempering simulation. Each temperature also has an associated, constant move amplitude $`\sigma _i`$ that is adjusted at the onset in order to have a reasonable acceptance ratio. We note that the swap move is very fast in computing time, since the current energy of each configuration is stored. The initial conditions are $`n`$ random arrangements of atoms. We do not thermalize the systems, since the parallel tempering swaps are very effective in arranging the configurations according to their energies. A typical energy trace for a parallel tempering run is shown in Fig. 9.

In general, one may be interested in averages of all the systems at the various temperatures, $`T_i`$. In this case, parallel tempering can be used to study phase diagrams, and the additional equilibration given by the swapping helps the systems at low temperature to sample the probability distribution effectively. We are specifically interested in the system with the lowest temperature, which will hop between likely zeolite structures until the one corresponding to the experimental sample is found. We monitor the structure at lowest energy, and we stop when all the atoms are 4-coordinated and the diffraction term indicates that a good match has been found.

The choice of temperatures is very important in parallel tempering. To determine them, it is useful to consider the energy fluctuations. We plot the energy histograms of the Monte Carlo data for a parallel tempering run, and we construct the temperature ladder so that the histograms overlap significantly. A typical example is shown in Fig. 10. A good choice for the temperatures can usually be obtained from an initial simulated annealing run. This allows us to locate the freezing temperature, a high temperature, and a low temperature. If necessary, the temperature selection can be refined to ensure that the energy histograms overlap and that the copies traverse the temperature ladder in the parallel tempering. This criterion leads to a ladder of temperatures that is automatically inferred from the properties of the system, rather than guessed or estimated. From Fig. 10, for example, it is clear that the high temperatures can be spaced more widely than can the low temperatures.

## IV Results

To assess the usefulness of this structure solution method we test it on the 118 publicly known zeolite structure types. For each of these materials, the chemical composition and the atomic positions are known. We describe here how the method fares on these known materials.

Since we do not have actual experimental diffraction information for most of these materials, we use the available data to construct synthetic X-ray diffraction patterns. The available data include information about the unit cell size and parameters, the space group symmetry, the type of atoms present in the unit cell, the atomic positions, the occupancies, and the Debye-Waller factors. The reflections included in our diffraction pattern are the ones in the range $`5^\mathrm{o}2\theta <35^\mathrm{o}`$. This excludes the low angle data, which usually have background contributions, and the high angle data, which are usually poorly resolved. The multiplicities of the peaks are accounted for in the production of the powder patterns. We use a wavelength of $`\lambda =1.54056`$ Å. We assume a crystallite size of $`1\mu `$m, which is typical for new zeolite samples. We assume that the data could be collected on synchrotron, and so we use instrument parameters appropriate for a beam line. With these assumptions, the diffraction pattern has a peak resolution of approximately $`\mathrm{\Delta }(2\theta )=0.06`$. We use this relatively conservative criterion when forming the composite peaks. Of course, the multiplicities of the peaks were accounted for in the production of the powder patterns. We verified that we could reproduce the intensities obtained with Cerius2, such as the ones in Fig. 2. In the cases where the framework topology can be described by a higher symmetry, we often, but not always, use the higher symmetry setting for the solution. This reduces the number of degrees of freedom, making the solution simpler and faster. As before, $`n_{\mathrm{unique}}`$ is the number of crystallographically distinct T-atoms used in the simulated annealing, and $`n_{\mathrm{max}}`$ is the number of crystallographically distinct T-atoms in the maximal symmetry setting. In the cases where several zeolite samples are available for a given framework, we generate the diffraction pattern using the material with the highest silicon content, or the lowest number of non-framework species.

The result of applying the solution procedure to the known zeolites is shown in Tables II and III. For each framework, we attempted a simulated annealing run with biased Monte Carlo moves. At the end of the run, we computed the coordination sequence of the unique atoms. The coordination sequence is a list of integers that counts the number of neighbors one, two, and so on connections away. It uniquely identifies a given structure through its topology, rather than through the precise locations of the atoms.

In the cases where we were unable to solve the structure in a few simulated annealing runs, we turned to parallel tempering. Using the energy histograms collected at the various temperatures in the simulated annealing run, we set up a ladder of five or six temperatures, using the corresponding move amplitudes, $`\sigma _i`$, from the simulated annealing. Using parallel tempering we were able to solve all the structures that we attempted. Depending upon the complexity of the structure, which is roughly proportional to the number of atoms in the unit cell, a solution is achieved in 0.2–4 hours on a Silicon Graphics Indigo<sup>2</sup> with a 195 MHz R10000 processor.

## V Discussion

From the results shown in Tables II and III we see that the introduction of biased moves in the simulated annealing dramatically improve the success rate of the method. The $`N_{\mathrm{MC}}`$ column refers to the number of simulated annealing attempts needed to solve the structure with simple Metropolis Monte Carlo moves, while $`N_{\mathrm{BMC}}`$ refers to the number of attempts required with biased moves. It is apparent that in most cases the biased moves substantially improve the sampling, allowing one to find the correct structure in fewer trials. More importantly, the technical limitations encountered with Metropolis Monte Carlo have been mostly removed, since most structures not solvable with simple moves can be solved with biased moves. It must be noted that several of the structures shown in the Tables were not known at the time Ref.  was published, and, arguably, some of the simpler new structures could have been solved with Metropolis moves. Indeed, at least four groups have used the Metropolis Monte Carlo approach, implemented in Cerius2, to solve new zeolite structures. Zeolite beta provides an example of a well known and important zeolite that can be solved with the biased Monte Carlo, but not with Metropolis Monte Carlo. The structures that were not solved with simulated annealing, and only these, were attempted with parallel tempering. All were solved. From the success of this approach, it is clear that parallel tempering is a much more powerful method than simulated annealing. We believe that parallel tempering is a much better method precisely because it samples the correct equilibrium distribution.

The most complex of the publicly known zeolites is ZSM-5. The framework structure type is MFI, and it is shown in Fig. 1. This structure is the most complex because it has the highest number of unique T-atoms, $`n_{\mathrm{unique}}=12`$. We were able to solve this structure in one day’s work that included the simulated annealing run, selecting the parallel tempering temperatures, and performing the actual parallel tempering run. The energy traces for the latter are shown in Fig. 9, and one can clearly see how the correct structure is found in one of the high temperature systems and then swapped down to the lowest temperature system. It is also clear from the trace that without the swapping, the lowest temperature system would never cross the energy barriers separating the initial condition from the correct structure. To showcase the power of the method, we solved solved NU-87, framework code NES, in a low symmetry setting with 17 unique atoms. This proves that parallel tempering is powerful enough that the use of the maximal symmetry setting is not required. Although an upper limit to the practical applicability of this method must exist, that limit is not obvious from the results of our trial runs.

We encountered a few types of problems when testing our method on the known zeolites. In some cases the framework has loops of length 3, with 3 T-atoms connected in a triangle. Clearly the bond angles in this case are far from the usual tetrahedral value, and the potentials that we defined in Section II may be incorrect for these particular structures. In other words, the correct structure may not correspond to minimum of the geometric terms in the figure of merit. In these cases, we found that a small adjustment of the geometric weights, $`\alpha _i`$ in Eq. (3), is a sufficient remedy. These cases can be recognized by a visual inspection of the structures produced, even with no a priori knowledge of the correct structure. The case of open frameworks, frameworks for which not all of the T-atoms are 4-coordinated, can be treated by not penalizing 3-coordinated atoms. These frameworks are listed in Table IV.

A few structures were harder to solve because of their merging. This is the case, for example, for MEP, PAU, and MWW. Let us consider MWW, which has $`n_{\mathrm{unique}}=8`$, $`n_{\mathrm{symm}}=24`$, and $`n_\mathrm{o}=72`$. The difficulty in this case stems from the fact that two combinations of merged atoms give the same $`n_\mathrm{T}`$: ($`12\times 5`$ atoms + $`4\times 3`$ atoms) $`=72`$ atoms and ($`12\times 4`$ atoms + $`6\times 4`$ atoms) $`=72`$ atoms. Many more combinations give numbers close to 72. The density term, $`H_\mathrm{D}`$, of Eq. (3) will not distinguish between these combinations. This occurrence can be inferred by visually inspecting the configurations of minimum energy generated or by tabulating the number of merged atoms associated with each unique atom. Again, no knowledge of the correct structure is needed to realize that the figure of merit has two or more deep wells separated by high barriers. The case of MWW is further complicated by the presence of 3-loops and the presence of two atoms which, although not connected, are separated by just 3.6 Å. In the case of MWW only, we adjusted the parameters $`\alpha _i`$ slightly, making $`\alpha _\mathrm{D}`$ small and reducing the range of the interactions. This allowed us to solve the structure with 3 iterations of parallel tempering. Because of these geometric irregularities, MWW was actually the most difficult structure for us to solve, even though it has only 8 unique T-atoms.

The presence of template molecules or heavy cations in the structure may cause a greater problem. In this case, the $`H_{\mathrm{PXD}}`$ term favors the presence of a nonzero scattering density in the regions occupied by the non-framework species. This makes the diffraction term ambiguous, and one can find many incorrect structures that are still feasible from the point of view of geometry alone. The preferred solution in this case is to calcine the structure to remove the template or to exchange the heavy cations for lighter ones. If the material is not stable to this treatment, then the non-framework species can be included as degrees of freedom in the figure of merit.

Zeolites with metal substitutions are quite common. Several of the structures that we solved are aluminophosphates or contain gallium, beryllium, cobalt, or zinc along with phosphorus. The effect of substitutions is quite dramatic, both on the geometry and on the diffraction. In fact, preferred bond lengths may change. The fact that the atomic species are different usually lowers the maximal symmetry allowed, and hence increases the number of unique atoms. We solved all known instances of these materials using just silicon to match the diffraction data. When possible, we used the higher symmetry allowed by assuming that all T-atoms are identical in the structure determination. Of course, it would be straightforward to simulate frameworks with different T-atoms and to include an atom exchange move in the Monte Carlo. This move, selected at random, would be accepted with a Metropolis criterion, and would lead to structures with the correct framework and the correct atomic species in each position. The case of the framework WEI, weinebeneite, is separate since this is a beryllium phosphate, and the bond lengths are quite different from regular aluminosilicates. A simple redefinition of the geometric terms in Fig. 4 would allow us easily to solve this material. A general extension would be to allow for species with with different coordination numbers as well as different T–T distance and T–T–T angle potentials. This extension would broaden the range of applicability of the method well beyond zeolites.

The issues of thermalization, autocorrelation times, and efficiency are usually addressed in numerical simulation studies. In our case, we are concerned with these elements in a qualitative sense only. In fact, all that matters for the utility of our method is that it determine the structure of an unknown material in a reasonable amount of time. We have shown that this time is very reasonable and very small compared to the time it takes to synthesize a new zeolite. The time taken by this approach should be compared to the time it takes to solve a single structure with conventional methods, typically measured in weeks, months, or year, or years. Timing issues may become relevant when one is interested in exhaustively listing all of the possible low energy geometries for a given cell size and symmetry, without the aid of experimental data. Hypothetical structures of this sort could help, for example, in the design of new zeolite materials, when the synthesis process is better understood.

## VI Conclusion

We have introduced a powerful biased Monte Carlo approach that can determine the structure of a new zeolitic material from the powder diffraction pattern and density, both of which are easily measured in experiments. All of the publicly known zeolites were solved in a realistic test application of the method. The method is rapid and automatic, making it a natural tool for use within the combinational chemistry paradigm. The proposed technique can also be used to generate hypothetical zeolite structures. The number of structures that can be generated in this way far exceeds the number constructed to date by hand. A database of structures can be built that will allow synthetic chemists to search for a structure that matches the X-ray powder diffraction pattern of newly made materials. The generation capability may also be directed, by, for example, specifying large pores, low-densities, or other interesting functional properties. Such use allows the creation of hypothetical zeolites with new, tailored structures. These structures can then be sought in rational syntheses. The direct, real space approach may also be applied to small, low quality molecular crystals, such as those from drugs, dyes, pigments, and organic non-linear optical materials, for which only powder diffraction data are available.

Perhaps of more general interest to the simulation community is the value of our work as a case study on the effectiveness of parallel tempering. We have shown that parallel tempering, combined with biased Monte Carlo, is a powerful method for molecular systems with many local minima. We have suggested a general and simple histogram method to determine the temperatures required in the extended ensemble. By considering how the energy histograms overlap, we generally expect that, away from critical points, the number of required temperatures is proportional to $`[E_{\mathrm{mix}}(N)E_0(N)]/\sqrt{N}`$, where $`N`$ is the number of degrees of freedom, $`E_0`$ is the energy at the temperature of interest, and $`E_{\mathrm{mix}}`$ is the energy at which the system can overcome all relevant barriers. How the number of required temperatures scales with system size will depend entirely on how $`E_{\mathrm{mix}}`$ scales with $`N`$. In simple cases, $`E_{\mathrm{mix}}`$ may scale as a correlation volume. In the worst case, $`E_{\mathrm{mix}}`$ will scale as $`N`$, and so the number of temperatures will be proportional to $`\sqrt{N}`$. By formulating the method as a Markov process on the level of a single move, we were able to provide a general and more efficient strategy for the choice of the individual system update frequencies. Based upon our experience with zeolite structure solution, we recommend use of parallel tempering whenever simulated annealing struggles on a minimization problem, especially if correct sampling of a low-temperature distribution is desired.

###### Acknowledgements.

It is a pleasure to acknowledge useful discussions with John Newsam and Giorgio Parisi. MF acknowledges useful discussions with Gerard Jungman. This research was supported by the National Science Foundation through grant CTS–9702403, by the University of California Energy Institute, by the Donors of the Petroleum Research Fund, and by Molecular Simulations, Inc.