## research papers

## Reverse Monte Carlo with geometric analysis – RMC+GA

^{a}Royal Institution, Albemarle Street, London, UK, and ^{b}Department of Earth Sciences, Downing Street, Cambridge, UK^{*}Correspondence e-mail: stephen@ri.ac.uk

An analysis of simulated framework structures based on the rigid-unit model quantifies the distortion of polyhedra from their ideal geometric forms and decomposes the motions of the structure into components of rigid-unit displacement, rigid-unit rotation and distortion. Case studies analysing the behaviour of quartz and of other silicates demonstrate that the method provides a novel way of extracting information from reverse Monte Carlo simulations. A program called *GASP* has been developed to perform this analysis and is freely available to researchers. Rotations are handled using the rotor method of geometric algebra.

Keywords: simulated framework structures; flexible framework structures; rigid-unit model; reverse Monte Carlo simulations; geometric analysis; computer programs.

### 1. Introduction

While crystal structures might be expected to be highly rigid, framework structures such as aluminosilicates in fact display considerable internal flexibility *via* cooperative motions of polyhedral groups of atoms. This fact has important consequences for the modelling of framework structures. We present a methodology, namely reverse Monte Carlo with geometric analysis, that quantifies the significance of such flexibility in the dynamics of framework structures.

The question of the rigidity or flexibility of framework structures in three dimensions remains poorly understood. Constraint counting as used by Maxwell fails to capture the flexibility of structures in which symmetry makes certain constraints redundant. This shortcoming is particularly relevant to crystal framework structures such as silicates, which possess many _{4} tetrahedra rotate and translate without distortion (Dove *et al.*, 2000). These rigid-unit modes (RUMs) can manifest themselves as soft modes in displacive phase transitions (Giddy *et al.*, 1993; Hammonds *et al.*, 1996), as sources of negative (Heine *et al.*, 1999; Welche *et al.*, 1998) and as mechanisms for large-amplitude thermal motion of vertex atoms at low energy cost (Gambhir *et al.*, 1999).

RUMs can be identified in *et al.*, 1994). While this approach can successfully identify candidate soft modes for phase transitions, and quantify to a degree the rigid-unit flexibility of a structure, the model does not capture the real-space motion of the structure. Firstly, this approach operates in the limit of infinitesimal amplitude; secondly, rotations are inherently anharmonic, as the motion of atoms is not linear in the rotation angle; and thirdly, rotations are not commutative. A motion of atoms due to a large number of active RUMs will therefore differ from a superposition of RUM eigenvectors.

Reverse Monte Carlo (RMC) modelling allows us to simulate configurations of atoms based directly on experimental data (McGreevy & Pusztai, 1988; McGreevy, 1995). The RMC technique used to generate the configurations in this paper was modified to account explicitly for both long-range and short-range order (Tucker *et al.*, 2000; Tucker, Dove & Keen, 2001; Tucker, Keen & Dove, 2001). Our configurations are consistent with the total scattering intensity, the pair distribution factors and the intensities of the Bragg peaks.

We have been working for some time on a method to extract information on RUM motion in real space from such techniques (Wells, Dove & Tucker, 2002; Wells, Dove, Tucker & Trachenko, 2002; Wells, 2003) in order to complement the reciprocal-space treatment of rigidity described by Hammonds *et al.* (1994). The principal problem is to model the behaviour of the polyhedra formed by real atoms in such a way as to relate it to the behaviour of the ideal polyhedra that are the objects of the rigid-unit model. Our analysis now not only quantifies the total atomic motion but also allows for both rotational and displacive motion of the polyhedra, whereas the earlier approach of Wells, Dove & Tucker (2002) and Wells, Dove, Tucker & Trachenko (2002) was restricted to rotational motion only.

We use the language of geometric algebra (GA; Lasenby *et al.*, 2000), a form of Clifford algebra, in which the basis vectors of three-dimensional space and Hamilton's quaternions exist within the same algebra. It is thus particularly well suited to the description of rotations. None of our results depends in detail on the form of the rotation operator used (it is simply computationally convenient for us to use the GA formalism) and so we will only quote those expressions necessary for the interpretation of our results.

The combination of RMC modelling with geometric analysis represents a new methodology, RMC+GA, for obtaining information on the behaviour of framework structures. This geometric analysis can also be applied to configurations generated by other methods, such as

(MD).### 2. Finding best-fit polyhedral motions

Suppose that we have two slightly different forms of the same structure, in which each polyhedron has undergone some displacement, rotation and distortion due to the motion of the atoms but there has been no change in topology, *i.e.* no bonds have been broken or have formed. Specifically, let us take two independently generated RMC fits to the same data; since RMC includes a stochastic element, these fits will differ in detail as to the positions of the atoms. As argued by Wells, Dove, Tucker & Trachenko (2002), we take each fit as an instantaneous snapshot of the system in dynamic disorder. Two independent (uncorrelated) fits can be taken as two snapshots separated by an infinite time interval.

We could also compare two different frames from a molecular dynamics simulation, either successive frames separated by a single time step or two frames with a wider temporal separation.

By comparison of the two structures, polyhedron by polyhedron, we can identify the displacements and rotations of the polyhedra, and the remaining distortion, namely the difference in shape between the two forms of the same polyhedron. We proceed by defining `mismatch scores', representing the differences between the two structures, and then minimizing these scores as a function of the rigid-unit motion of the polyhedra. We can also compare the polyhedra in a simulated structure with those in an idealized structure with perfectly regular polyhedra, so as to quantify the distortion of the polyhedra or to identify an RUM acting as a soft mode for a phase transition.

The steps in the comparison process are as follows: comparison of atomic positions between the two configurations; displacement of polyhedra so that their centres coincide; and rotation of the polyhedra so that their orientations match as closely as possible. The process is illustrated in Fig. 1. Effectively, we decompose the motions of all the atoms of the polyhedron into a bodily displacement (given by the displacement of the central atom), a bodily rotation and a residual distortion.

#### 2.1. Definition of mismatch scores

We define an initial mismatch score to be the sum of the squares of the displacements of all the atoms from one configuration to the other. Thus, if the atoms in one configuration have positions *r*_{i} and each atom's position in the second configuration is , the initial mismatch score before any fitting takes place is now defined to be

Since this score is simply a sum over atomic contributions, we can either consider the mismatch score for the structure as a whole or decompose the score by atomic species; for each species, we obtain a mean-square atomic motion

and an r.m.s. atomic displacement

while, of course, the r.m.s. atomic displacement taken over all atoms is simply

Since *M*_{before} scales linearly with the size of the system, comparison between different systems is easier if we divide *M* by the size of the system – that is, the number of polyhedra – to obtain the value of *M* per polyhedron. The mismatch values in the case studies that follow are given in this form.

Having obtained these measures of the atomic motion, we now consider the structure as an assemblage of polyhedra, decomposing it as illustrated in Fig. 1. We make the approximation that the displacement of each polyhedron is equal to the displacement of its central atom. We now wish to identify the rotational motions of structural polyhedra. Let us consider the structure bond by bond. For each bond, there is a vector from the central atom to an atom *q* at a vertex of the polyhedron *p*; for one form of the structure, we call this vector **pq**, for the other form, . Our assumption, from the rigid-unit picture, is that there exists a rotation of the polyhedron *p* that takes the set of vectors very close to the set . We can write down a vector , representing the vector **pq** after a rotation *B*; this allows us to write a mismatch vector, **m**, given by . Minimizing the total magnitude of all the mismatches for the polyhedron, , with respect to the parameters of the rotation *B* will give the rotation that best fits **pq** onto . The value of *M*_{p} before minimization with respect to *B* represents the degree of mismatch between the polyhedra when they have been displaced so that their centres coincide but they have not yet been allowed to rotate. We therefore designate this value as *M*_{p,fit centres}. The final value of *M*_{p}, representing the residual distortion of the polyhedra, we write as *M*_{p,after}.

*M*_{after}, the average value of *M*_{p,after} for all the polyhedra in the structure, can be compared with the original measure of the degree of atomic motion, *M*_{before} per polyhedron, in order to obtain a measure of the importance of rigid-unit motion. For concision, we frequently refer to *M*_{before} as the disorder score and to *M*_{after} as the distortion score.

Full mathematical details of the process of finding the best-fit rotation are given by Wells, Dove & Tucker (2002). However, note that the quantity defined in that paper as `*M*_{before}' or as `total disorder' is, in fact, the sum over all polyhedra of our *M*_{p,fit centres} and is not related to *M*_{before} or the disorder score as defined in the present paper.

#### 2.2. Interpretation of rotors

We describe a rotation using three independent variables, *B*_{x,y,z}, which define a `bivector' object describing a `rotor' operator; these are concepts of geometric algebra. To interpret the results in the case studies that follow, it suffices to interpret *B* in terms of an axis and an angle . The size of the rotation angle , in radians, is given by

The axis of the rotation is given by a unit vector, , whose Cartesian components, *b*_{i}, are given by

The case *B* = 0 is the identity operation, that is, *R*(0) = 1. To first order, therefore, the values of *B*_{x,y,z} are the rotations about the *x*,*y*,*z* axes, in radians.

#### 2.3. How important is rigid-unit motion? Comparison of mismatch scores

We wish to define a measure of the importance of rigid-unit motion in the dynamics of a structure. Perfect accommodation (zero distortion score *M*_{after}) cannot be this criterion, since even when RUMs are dominant, higher-frequency distortive modes will also be active. Some distortion of the polyhedra is therefore inevitable. In the disorder score, each atom contributes the square of its displacement, . For each atomic species we find an r.m.s. atomic displacement

The distortion score is the sum of the mean-square displacements of the vertex atoms from their geometrically ideal positions; when different types of polyhedra are present, this value will, of course, be different for each. Thus, from the distortion score, we can obtain the r.m.s. distance of vertex atoms from ideal vertices, *d*_{atom vertex}, and these distances can be compared with the r.m.s. atomic displacements in order to gauge the importance of rigid-unit motion. Our identification of the displacement of central atoms with the displacement of polyhedra means that central atom species make no direct contribution to the distortion score. The significant comparison is therefore between the motion of vertex atom species and the distortion of polyhedra. In the case of quartz discussed below, for example, the r.m.s. motion of O atoms is twice as large as the r.m.s. distortion, *d*, at room temperature, and four times as large at high temperatures in the phase.

The value of *M*_{fit centres} is not particularly significant in itself. We do, however, note that, when rotational amplitudes are large, *M*_{fit centres} lies at or above the value of *M*_{before}. If, on the other hand, displacements are more significant than rotations, *M*_{fit centres} will lie below *M*_{before}. Thus far, we find rotational-dominant behaviour in the crystalline framework silicates such as quartz, while silica glass also shows significant displacive motion.

Evidently, the fitting process generates a rotor for every polyhedron in the structure. The r.m.s. magnitude of the rotor, or of a component of the rotor, provides an indication of the magnitude of rigid-unit motion in the structure, and it may be of interest to observe this magnitude as a function of temperature. Furthermore, it is possible to plot the distribution of rotor components for all the polyhedra, rather than finding an r.m.s. value; this process may yield information on the nature of the rigid-unit motion. We use both forms of analysis in our discussion of quartz below.

### 3. Applying a force model *via* geometric modelling

The geometric model on which this analysis is based, in which each polyhedron exists both as a set of real atoms and as an ideal geometric form, can itself be used to generate configurations of atoms by establishing the ideal forms first and then constraining the real atoms to fit them. This process provides a simple means of studying the geometric response of frameworks to substitutional defects or compression, and of rapidly relaxing idealized structures (for example, silica networks with 180° Si—O—Si angles) to obtain realistic bond lengths and angles.

The model requires several iterations over two distinct steps. In the first step, we fit a geometrically ideal polyhedron over each real polyhedron using the GA rotor-fitting approach, so that the vertices of the ideal polyhedra provide a close match to the positions of the vertex atoms. The initial fitting of the ideal polyhedron involves only three variables (the parameters describing the orientation of the polyhedron), and this step is easily performed using the method given by Wells, Dove & Tucker (2002). In the second step, each atom moves within a simple-harmonic force model based on the positions of the vertices. The displacement of the atom relative to the position of an ideal vertex can be divided into the components parallel and perpendicular to the direction of the ideal vertex–centre (*e.g.* Si—O) bond. The parallel component represents the extension of the bond, while the perpendicular component represents distortion of the internal bond angles of the polyhedron. To each component we associate a spring force, tending to pull the atom towards the ideal position. The three-body terms constraining the O—Si—O bond angles are thus replaced by constraints on the displacement of the vertex atoms from their ideal positions. In keeping with the principle of using two-body rather than three-body potentials, the centre–vertex–centre (*e.g.* Si—O—Si) bond angle constraints are represented as constraints on the centre–centre distances.

Having relaxed the atoms within this simple force model, we then refit the ideal polyhedra over the real polyhedra using rotor fitting and iterate until the structure is fully geometrically relaxed. The geometrically ideal polyhedra therefore represent a form of multibody potential in which the evaluation of bond angles is unnecessary, as angular information is implicit in the shape of the polyhedron. It is, of course, possible to include constraints on the relaxation by, for example, pinning certain atoms or polyhedra in place. This allows for a direct connection between the real-space and reciprocal-space pictures of polyhedral motion, since the ideal polyhedra of the model are the objects of RUM theory. Thus we can both obtain the RUM component of a given set of atomic motions from the motions of the ideal polyhedra, and apply one or more RUM eigenvectors to the ideal polyhedra and observe the response of the real polyhedra.

This geometric modelling process has been used successfully to model the high-*P* deformation mechanism of a zeolite framework with EDI topology (Gatta & Wells, 2004).

### 4. Program availability and capability

We hope that the analysis will be of interest to the crystallography community. We have therefore implemented it in a self-contained program, called *GASP* (geometric analysis of structural polyhedra). The source code is written entirely in Fortran 90 and is available with documentation from http://www.esc.cam.ac.uk/minsci/downloads/gaspfiles.tar.gz . Questions, comments, bug reports and fixes should be reported to the corresponding author. The program in its current version (1.5) performs three principal functions. The first function, known as preconditioning, analyses a structure and identifies the polyhedra within it. During this process a connectivity file is generated, recording which atoms are bonded together, and an ideal polyhedron is fitted over each real polyhedron, quantifying the distortion of the polyhedra from their geometrically ideal forms. If the connectivity is already known, it can be provided as input. Tetrahedra, octahedra and planar triangles can all be present.

The comparison function takes as its input two structures that have been preconditioned, that is, decomposed into their component polyhedra. This function performs the analysis described in §2.1, quantifying the difference between the two structures and analysing it into components of rigid-unit displacement, rigid-unit rotation and distortion. The two structures could be two independent RMC fits, an ideal (average) structure and a fitted (disordered) structure, two frames from an MD simulation, or representations of a structure before and after the introduction of a defect.

The third function makes use of the modelling technique described in §3 to relax a structure towards geometric ideality. The program and analysis are designed for the comparison of structures containing very large numbers (of order 10^{2}–10^{5}) of atoms. The program is driven by a text file containing geometric information and program commands. Structures can be given either in a native .dat format or in the common .xtl format. Full information on file formats and commands is available with the code.

### 5. Case study

Our principal case study is the analysis of a series of quartz configurations generated by reverse Monte Carlo fitting to total neutron scattering data. Frameworks of corner-linked SiO_{4} tetrahedra such as these are the classic objects of the RUM model, and such frameworks are therefore a suitable first application for our real-space analysis. Our analysis provides information on the RUM contribution to the dynamic disorder and the effects of phase transitions on the RUM behaviour. We quantify the significance of RUM motions and show that they are very important in the dynamic disorder, particularly at high temperatures. By analysis of the and phases, we show that the disorder in the high phase can be explained on the basis of the excitation of RUMs. We show the effect of the on the RUM behaviour in the phase, demonstrating the link between the symmetry change due to the transition and the amplitudes of rotations about different axes.

While quartz shows considerable geometric flexibility, it is not the most flexible crystalline silicate. We compare our results for quartz with those for cristobalite and zeolite Y, demonstrating that larger amplitudes of atomic motion are obtained in the more RUM-flexible frameworks.

#### 5.1. Quartz

In a recent total neutron scattering study of quartz (Tucker *et al.*, 2000; Tucker, Keen & Dove, 2001), measurements were performed over a wide range of temperatures below the temperature and for some temperatures above the The RMC configurations therefore provide information about changes in both short-range and long-range order as a function of temperature on heating through the The configurations contained 6000 SiO_{4} polyhedra. Sections of some of the configurations so generated are shown in Fig. 2. It is clear that, as temperature increases, the variation about the `average' structure (the structure found from the Bragg peaks) increases markedly. Since RUMs provide a mechanism for large reorientations of the polyhedra at small energy cost, we would expect RUMs to play an important role in this variation. By visual inspection, it would seem that the polyhedra are undergoing rigid-unit motion. The rotational motion and polyhedron distortion were considered by Wells, Dove, Tucker & Trachenko (2002). We can now complete that analysis by considering also the total atomic motion. Fortunately, data were available for a large number of points covering both the low and the high phase, and Fig. 3 shows the degree of disorder in quartz as a function of temperature. This graph shows the mismatch scores defined in §2, in Å^{2} per polyhedron, between two independent RMC fits to the scattering data at each temperature.

Note that the tetrahedral distortion is much less than the total disorder. The large difference between the two indicates the contribution of RUMs to the dynamic disorder. The value of *M*_{fit centres} lies slightly above the value of the distortion score *M*_{before}, indicating that the RUM motions are predominantly rotational rather then displacive, as discussed in §2.3. Note particularly that the distortion score is relatively insensitive to temperature and to the while the total disorder increases sharply with temperature and shows the quite clearly.

The tetrahedral distortion can be further divided into the terms due to variation in bond lengths and due to bending of the O—Si—O bond angle. It is clear that the distortion is primarily due to bond bending, with a much smaller bond-stretching component. The nature of the distortions does not appear to change over the temperature range of this study. Wells, Dove, Tucker & Trachenko (2002) showed that there is no correlation between the distortions in the two different configurations.

Fig. 4 shows the breakdown of the mismatch scores into r.m.s. motion values for Si and O atoms, plus distortion values for vertex O atoms in SiO_{4} polyhedra. The motion of O atoms is twice the distortion value at room temperature and in the phase reaches four times the distortion value. This comparison is very revealing of the importance of rigid-unit motion. Since the displacement of the Si atoms is identified with the displacement of the polyhedra, they do not contribute directly to the distortion value.

#### 5.2. in quartz

In terms of the average structure, the phase of quartz is obtained from the phase by the condensation of an RUM as a soft mode (Vallade *et al.*, 1992). This RUM involves rotations of the polyhedra about the axes, with the rotation angle reaching ∼17° at room temperature. It should be possible to detect this rotation using geometric analysis. We therefore compare the fitted configurations with an ideal (zero-order parameter) high-temperature form and observe the distribution of rotations of the fitted polyhedra. We expect to observe rotations of the polyhedra about the axes, which in terms of GA means the following; if we plot the *X* and *Y* components of the rotors for the polyhedra, relative to the high phase, against each other, the results should cluster about three points. One point will be at (where radians or 17° at room temperature), while the other two points will lie at 120° to the first point in the *XY* plane. If more than one domain were present in the sample, there would be six equally spaced clusters instead of three.

This prediction is amply born out by the results given in Figs. 5–8.

With increasing temperature, the distribution of rotations relative to the high phase becomes broader as the rigid-unit component of the dynamic disorder increases in amplitude. In the phase, the distribution of rotations becomes isotropic at about zero and is very broad. It appears that the disorder in -quartz arises largely from the excitation of the RUM spectrum.

The value of the tetrahedral rotation angle, , as obtained by fitting Gaussians over the rotor distributions and extracting the peak position, is plotted in Fig. 9. At room temperature, we find an average rotation of 17° relative to the phase, as expected.

| Figure 9 Unfortunately the distribution of angles near the transition is so broad that the parameter cannot be obtained with any accuracy in this region. |

##### 5.2.1. Amplitude of rotational motion

The – . For each polyhedron, there is a rotation describing its reorientation from one configuration to the other. The distribution of rotor components *B*_{x,y,z} describing these rotations is a Gaussian centred on zero. The r.m.s. value of the rotor component is a measure of the amplitude of rotational motion, as shown in Fig. 10. It appears from the behaviour of the rotor components that, in the phase, the polyhedra are freer to rotate about axes lying in the *xy* plane than they are about the *z* axis, a fact that reflects the hexagonal symmetry of the structure. The effect of the however, with its imposition of an average rotation of the polyhedra about the axes, is to constrain the *x* and *y* rotations and reduce them to the same amplitude as rotations about the *z* axis.

##### 5.2.2. RUM frequencies and amplitudes and their relation to the order parameter

We now look more closely at this effect, whereby the symmetry-breaking in the

affects the rotational behaviour, by considering the link between the the frequencies of the RUMs and their amplitudes. A simple Einstein model for the rotational oscillations of the polyhedra suggests thatand we therefore plot , a measure of , against *T* to observe the behaviour of the RUM frequencies. The split-atom model (Giddy *et al.*, 1993, Hammonds *et al.*, 1994) would suggest that, as the average tilt, , of the polyhedra relative to their positions in the phase increases from zero, modes that are RUMs in the phase gain frequency in proportion to . Since the in quartz is of almost tricritical character (Carpenter *et al.*, 1998), and the data do not approach the closely enough to allow us to observe any first-order discontinuity, we take it that

As a simple model, we suppose that the value of for the high phase is a result of a population of RUMs having a typical frequency , and that . We extrapolate linearly from the values in the high phase to obtain the quantity , that is, the predicted value of if all the RUMs from the high phase remained RUMs in the low phase.

We now suppose that RUMs in the low phase develop a component of tetrahedral distortion and gain frequency in proportion to the growth of the

. These modes will therefore have their amplitude reduced and our prediction for then becomesWe perform two separate fits, one to the average of the *X* and *Y* rotor components, and another to the *Z* rotors. The results are plotted in Fig. 11. For the *X* and *Y* rotors, the model fits the data well, with a value of *K* = 0.75. For the *Z* rotors, however, the fit is best with *K* = 0.19. This outcome is understandable in terms of our model in which rotations about the *z* axis are more constrained than those about the *x* and *y* axes in the high phase.

This result confirms the impression obtained from the rotor components in Fig. 10; rotations about the *z* axis are constrained throughout, while rotations about axes in the *xy* plane are freer in the high phase but become constrained because of the symmetry change in the phase transition.

#### 5.3. Comparison with other crystalline silicates

As we have seen, quartz displays a considerable degree of geometrical flexibility; however, it is not the most flexible of the crystalline silicates.

Tucker, Squires *et al.* (2001) demonstrated that the structural disorder in the high phase of cristobalite could be accounted for on the basis of rigid-unit motion, without the need for a domain model. Wells, Dove, Tucker & Trachenko (2002) considered the flexibility of the cristobalite structure from the point of view of RUM theory, showing that it is greater than that of quartz.

Zeolites are another class of framework structures for which an RUM analysis is fruitful. The high symmetry and large unit cells of some zeolite structures give them a large number of RUMs, and in some cases, the number of RUMs can even be a finite fraction of the total number of modes (Hammonds *et al.*, 1998).

We have performed geometric analyses on cristobalite configurations fitted to total neutron scattering data at five temperature points and on siliceous zeolite Y configurations fitted to total neutron scattering data at two temperature points, using the same RMC technique as for the quartz. A geometric analysis of the cristobalite configurations, considering only the rotational motions, was performed by Wells, Dove, Tucker & Trachenko (2002). We now present the results of our full analysis in the form of the r.m.s. motions of the Si and O atoms, and the r.m.s. distortion of the SiO_{4} polyhedra (Fig. 12). The corresponding values for quartz, from Fig. 4, are also given for comparison.

The zeolite and the cristobalite display very similar levels of atomic motion and of tetrahedral distortion; these values may therefore be typical of highly flexible crystalline silica. The degree of atomic motion is considerably higher than that seen in quartz, reflecting the greater geometric flexibility of these two structures. For the zeolite, cristobalite and quartz in the phase, the ratios of the r.m.s. oxygen motion to the r.m.s. tetrahedral distortion are similar and indicate the dominance of rigid-unit motion in the dynamic disorder, while in quartz at lower temperatures, this ratio decreases, as the symmetry change from the

decreases the geometric flexibility.### 6. Conclusions

We have developed an analysis connecting the motion of atoms in real space to the concepts of rigid-unit theory. Using this approach, we can quantify the importance of rigid-unit motion in dynamic disorder. We have shown that the method of `reverse Monte Carlo with geometric analysis' (RMC+GA) is extremely useful for extracting information from RMC simulations of framework structures. The method is implemented in a Fortran 90 code, *GASP*, which is freely available.

### Acknowledgements

The authors thank Kostya Trachenko, Alexandra Simperler and Diego Gatta for helpful discussions and suggestions. Funding for the thesis upon which this work is largely based was provided by the EPSRC.

### References

Carpenter, M. A., Salje, E. K. H., Graeme-Barber, A., Wruck, B., Dove, M. T. & Knight, K. S. (1998). *Am. Mineral.* **83**, 2–22. Google Scholar

Dove, M. T., Trachenko, K. O., Tucker, M. G. & Keen, D. A. (2000). *Rev. Mineral. Geochem.* **39**, 1–33. CrossRef Google Scholar

Gambhir, M., Dove, M. T. & Heine, V. (1999). *Phys. Chem. Miner.* **26**, 484–495. Web of Science CrossRef CAS Google Scholar

Gatta, G. D. & Wells, S. A. (2004). *Phys. Chem. Miner.* Submitted. Google Scholar

Giddy, A. P., Dove, M. T., Pawley, G. S. & Heine, V. (1993). *Acta Cryst.* A**49**, 697–703. CrossRef CAS Web of Science IUCr Journals Google Scholar

Hammonds, K. D., Dove, M. T., Giddy, A. P. & Heine, V. (1994). *Am. Mineral.* **79**, 1207–1209. CAS Google Scholar

Hammonds, K. D., Dove, M. T., Giddy, A. P, Heine, V. & Winkler, B. (1996). *Am. Mineral.* **81**, 1057–1079. CAS Google Scholar

Hammonds, K. D., Heine, V. & Dove, M. T. (1998). *J. Phys. Chem. B*, **102**, 1759–1767. Web of Science CrossRef CAS Google Scholar

Heine, V., Welche, P. R. L. & Dove, M. T. (1999). *J. Am. Ceram. Soc.* **82**, 1793–1802. CrossRef CAS Google Scholar

Lasenby, J., Lasenby, A. N. & Doran, C. J. L. (2000). *Philos. Trans. R. Soc. London Ser. A*, **358**, 21–39. CrossRef Google Scholar

McGreevy, R. L. (1995). *Nucl. Instrum. Methods Phys. Res. Sect. A*, **354**, 1–16. CrossRef CAS Web of Science Google Scholar

McGreevy, R. L. & Pusztai, L. (1988). *Mol. Simul.* **1**, 359–367. Web of Science CrossRef Google Scholar

Tucker, M. G., Dove, M. T. & Keen, D. A. (2000). *J. Phys. Condens. Matter*, **12**, 723–730. Web of Science CrossRef Google Scholar

Tucker, M. G., Dove, M. T. & Keen, D. A. (2001). *J. Appl. Cryst.* **34**, 630–638. Web of Science CrossRef CAS IUCr Journals Google Scholar

Tucker, M. G., Keen, D. A. & Dove, M. T. (2001). *Mineral. Mag.* **65**, 489–507. Web of Science CrossRef CAS Google Scholar

Tucker, M. G., Squires, M. D., Dove, M. T. & Keen, D. A. (2001). *J. Phys. Condens. Matter*, **13**, 403–423. Web of Science CrossRef CAS Google Scholar

Vallade, M., Berge, B. & Dolino, G. (1992). *J. Phys.* I, **2**, 1481–1495. CrossRef CAS Google Scholar

Welche, P. R. L., Heine, V. & Dove, M. T. (1998). *Phys. Chem. Miner.* **26**, 63–77. Web of Science CrossRef CAS Google Scholar

Wells, S. A. (2003). PhD thesis, University of Cambridge, UK. Google Scholar

Wells, S. A., Dove, M. T. & Tucker, M. G. (2002). *J. Phys. Condens. Matter*, **14**, 4567–4584. Web of Science CrossRef CAS Google Scholar

Wells, S. A., Dove, M. T., Tucker, M. G. & Trachenko, K. (2002). *J. Phys. Condens. Matter*, **14**, 4645–4657. Web of Science CrossRef CAS Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.