Obtaining diffuse scattering patterns from computer simulations – a retrospective

The paper describes how the calculation of diffuse scattering from atomistic model crystals has developed over the last approximately 50 years. Not only has the quality of observed diffuse X-ray scattering data improved immensely with the advent of electronic area detectors and synchrotron radiation but the enormous increase in computer power has enabled patterns, of comparable quality to the observations, to be calculated from a Monte Carlo model.

1. Introduction

In view of Hans-Beat Bürgi's interest in molecular crystals, their atomic displacement parameters, different kinds of disorder and diffuse scattering and computer simulation [see for example, Bürgi & Dunitz (1970), Capelli et al. (2000), Birkedal et al. (2003), Burgi et al. (2005), Nemkevich et al. (2010) and Michels-Clark et al. (2013)], we have chosen as our contribution to this special issue an account of our own journey through this fascinating field, beginning in the 1970s when photographic film was used to record diffuse scattering through to the present day where measurements using area detectors and synchrotron radiation are quite routine.

It is now well accepted that computer simulation of a model crystal provides a general method by which diffuse scattering of all kinds and from all types of materials can be interpreted and analysed. What it has been possible to achieve at any point in time has depended strongly on the power of available computers. During the last 40 years (since about 1980) there has been an increase in computational power by a factor of approximately 10⁵. Fig. 1 shows a plot of the evolution of computer power per cost in MIPS/$1000 (1997 dollars) (Note: 1 MIP = 10⁶ instructions per second).

Figure 1 Moore's Law. The plot shows the evolution of computer power per cost in MIPS per $1000 (1997 dollars). Since 1980 there has been an increase by a factor of approximately 10^5.

It is interesting to note that in 1980 when the first work described here was published the IBM PC, the Commodore 64 and the Macintosh 128 had not yet emerged. By present day standards computing was in a primitive state. Having said that it is salutary to realise that the Voyager I & II spacecrafts were launched in 1977 and their computers are still functioning in outer space today.

Computers have played a major role at all stages of our journey through the field of diffuse scattering and disorder. Computers are involved in three different aspects of the field:

(1) Data acquisition, data processing and equipment control.

(2) Model building and simulation.

(3) Calculation of diffuse scattering patterns from a model.

In the sections that follow we chart the progress that has been made by showing examples of what has been achieved at various stages. We begin with a study published in 1980 (Welberry & Jones, 1980) and move towards our most quantitative study to date that was published in 2014. Finally we discuss the recently published ultrafast method for calculating diffuse scattering from an atomistic model (Paddison, 2019) which promises to revolutionise the way in which studies may be carried out in the future.

2. Use of optical transforms

The original idea of Sir Lawrence Bragg of exploiting the analogy between optical diffraction and X-ray diffraction in solving structural problems received its greatest stimulus from the work of the Manchester School of Research under Professor H. Lipson in the early 1950's (Chaudhuri, 1976). With their background of experience, optical transform methods subsequently evolved as a useful aid in the elucidation of structures from X-ray diffraction patterns (Harburn et al., 1975).

In this section a brief account is given of a study of disorder in 2,3-dichloro-6,7-dimethylanthracene (DCLDME) in which optical transforms were used (Welberry & Jones, 1980). Fig. 2 shows a comparison of observed and calculated diffraction patterns for the h0l section of DCLDME. The observed X-ray pattern was recorded on X-ray film and digitised using a film scanner/writer [an Optronics P-1700 Photomation system (Harburn et al., 1974)] using a grid size of 100 µm. The same device was also used for rewriting the final reciprocal-lattice image after removal of the Weissenberg distortion. This is shown in Fig. 2(b).

Figure 2 Comparison of the diffraction patterns for the b-axis projection of the crystal structure of 2,3-dichloro-6,7-dimethylanthracene: (a) average structure drawing, (b) X-ray pattern and (c) optical pattern.

The same Photomation system device was used to plot on film a model of the disordered crystal structure. The Photomation is basically used as an accurate small-scale plotter. Data points may be written on a 12.5 µm (= 1 raster) grid over the area of the film (maximum width is 120 mm). The chosen scale for the plots was 10 raster units to 1 Å and this provides the following specifications for the plotted representation of the crystal structure. Since one raster unit (12.5 µm) represents 0.1 Å, positional coordinates are plotted to an accuracy of ±0.05 Å. The total size of the plot that is conveniently fitted on the 120 mm film is about 9000 raster units representing a linear dimension of ∼900 Å or 100 9 Å cells. With ∼10 000 unit cells, reasonable representations of statistical parameters could be achieved. The optical diffraction pattern from the model crystal of DCLDME is shown in Fig. 2(c).

From this study it was possible to estimate the magnitude of the correlation between the orientation of neighbouring molecules within the close-packed molecular ab planes as well as between neighbouring planes along c.

Although the Photomation was a sophisticated instrument relying on an optical encoder and accurate stepper motors to achieve its precise raster grid it required only a modest mini-computer to control it. Model coordinates were obtained on a separate mini computer (mostly a PDP-11) and transferred via magnetic tape. The models themselves were produced using simple programs that used growth disorder models involving both occupancy and displacement variables. These stochastic models (as their name suggests) were grown sequentially and were very fast and computationally undemanding [see Welberry & Galbraith (1973) and Welberry et al. (1980)]. They did not involve any iterative process as, for example, did Monte Carlo simulation.

3. Diffuse scattering studies circa 1990–2010

By 1990 a number of things had changed: the replacement of film by electronic detectors for recording diffuse X-ray scattering; the development of more sophisticated ways of building models of the disordered crystal; and the move away from optical transforms to direct calculation of the model diffraction patterns.

3.1. Measuring diffuse X-ray patterns

The recent development of position-sensitive counters (1D linear or even 2D area detectors) offered the possibility of replacing the use of film for recording diffuse scattering and allowing routine high-quality measurements to be made in the ordinary laboratory, using conventional X-ray generating equipment. These new counters gave only a limited advantage for measuring Bragg intensities (since advantage is gained only when several peaks are simultaneously incident on the detector), but they were ideally suited to the measurement of continuous diffuse distributions where all regions of the detector could be usefully utilized the whole time. Fig. 3(a) shows the system built in our laboratory which utilizes the curved PSD manufactured by STOE (Wölfel, 1983). Fig. 3(b) shows an example of the high-quality diffuse scattering data that could be obtained with this instrument.

Figure 3 (a) The position-sensitive detector (PSD) system used to record X-ray data from about 1990. (b) An example diffraction pattern recorded on the PSD system. The 0kl section of Bemb2 (1,3-dibromo-2,5-diethyl-4,6-dimethylbenzene) recorded at 100 K.

The open geometry of this system allowed easy access for a cryostream to be used for low-temperature studies. The collimation and beamstop system were designed to bring the incident X-rays to within a few millimetres of the sample, thereby virtually eliminating problems with air-scattering. See Osborn & Welberry (1990) for further details.

3.2. Construction of model crystals

Computers continued to evolve rapidly allowing models to make use of more sophisticated methods of simulation. The growth disorder models used earlier were now largely replaced by Monte Carlo simulation which involved many cycles of iteration to arrive at a disordered distribution close to an equilibrium structure. The quality of the observed patterns was such that it became necessary to consider how large a model crystal was needed to allow the calculation of a pattern of comparable quality. Two main criteria needed to be satisfied. The size of the model crystal needed to be large enough so that the resolution in real space allowed sharp diffuse features in the observed patterns to be satisfactorily captured. Here the linear dimension of the model crystal was a factor. Any correlations in the structure of a length greater than half the crystal size could not be captured. In addition the variations of intensity due to statistical fluctuations of the disordered patterns needed to be sufficiently small. This is proportional to N^1/2, where N (= N_x × N_y × N_z) is the total number of unit cells in the structure. Three-dimensional models could now be used but computer memory sizes were still very small by modern standards. One megabyte of RAM was still a pipe dream. A crystal size of 32 × 32 × 32 unit cells was adopted as a rule-of-thumb guide to simulation size. Fig. 4 shows example calculations for a simple disordered crystal. It is seen that the calculation made with the full 32 × 32 × 32 unit cells in Fig. 4(d) shows smooth and clearly defined diffuse peaks while the calculations made with smaller sizes do not.

Figure 4 Different quality diffraction patterns calculated from a model of the Tl correlations in TlSbOGeO4 (Welberry & Mayo, 1998). (a) The three types of correlation vectors. (b) Pattern calculated from a single 10 × 10 × 10 sample taken from a larger crystal. (c) The best pattern obtainable from a crystal of 10 × 10 × 10 unit cells. (d) The best pattern obtainable from a crystal of 32 × 32 × 32 unit cells.

3.3. Calculation of diffraction patterns

Optical transforms continued to be used through the 1980s and it was not until about 1990 that it became feasible for the first time to compute the diffraction patterns directly from the model disordered structures. A computed diffraction pattern would offer a number of advantages over the optical method. For example, in the optical method atomic scattering factors could only be represented fairly crudely by the use of different-sized apertures, the accuracy of positioning of the atoms was also rather limited (by the resolution of the film-writing device) and furthermore intensities were only obtained on an arbitrary scale. But perhaps the most important limitation of the optical method was that it was strictly limited to two dimensions (2D), and all that could be obtained from a 3D simulation was a 2D projection of a relatively thin slice of the 3D structure. A further constraint in carrying out these calculations was that the sampling in reciprocal space needed to be of a comparable resolution to that obtained in the PSD X-ray measurements (see §3.1). It was estimate that the resolution was such that a typical 2D section of data contained ∼45 000 independent data points.

3.4. Development of the program Diffuse

The total complex scattered amplitude from a crystal can be written as a simple sum of plane waves:

where

is the structure factor of the mth unit cell. is the position vector of the mth cell, r_n is the location of the nth atom in the mth cell, f_n is the atomic scattering factor of the atom n, k is the diffraction wave vector, N_m is the number of atoms in cell m and N is the number of unit cells in the crystal.

There are two obvious approaches that can be used to obtain the diffuse intensity of a simulated crystal using digital computers: compute the sum in equation (1) directly for each reciprocal coordinate desired, or use fast Fourier transform (FFT) techniques. The FFT is simply an algorithm that can be used to calculate discrete Fourier transforms quickly. In order to use the FFT all atoms in the crystal must fall on a uniform grid. The unit cells define a natural grid but only atoms on the corners of the cells would actually lie directly on it. In order to include all atoms in the cell, the basic crystal lattice must be subdivided further to bring each atom site as close as possible to a grid point. If each unit cell is divided into a 32 × 32 grid then placement of the atoms will be accurate to about 3% of the unit-cell dimension. This level of accuracy is only comparable to that provided by the film-writing device used in the optical transform method.

Consideration was given to the possibility of using FFT but although some advantage in speed could theoretically be obtained, the method required far more memory than computers of that era could typically provide. Consequently a direct method of computation using equation (1) was adopted. This led to the development of the computer program called Diffuse (Butler & Welberry, 1992). Direct Fourier summation over the atoms in the simulated crystal proved to be the most general and straightforward method of computation. Factoring the atomic scattering factors from the summations and careful design of the computational procedure provided the extra speed necessary to make this method viable. For example, instead of computing the large number of complex exponentials that are involved throughout the calculation, the values were first computed and stored in memory and then subsequently retrieved from a look-up table when required. Even with all these efforts to speed up the calculation it was found that in 1992 transforms comparable to those obtained optically required approximately 10 to 50 h using a VAXstation 3100 computer. Tests on a vector processing Fujitsu VP100 supercomputer showed the code to be easily vectorized and a speed improvement over the VAXstation 3100 of a factor of 100 was obtained.

It was found that averaging the intensities from a number of small regions or, `LOTS', inside the larger simulation produced superior results to those obtained by making a single calculation from the whole simulation. It is only necessary that the size of the `LOTS' is larger than the extent of the correlations. Fig. 5 shows schematically how the `LOTS' may be chosen randomly from the whole simulation. The quality of the calculated pattern increases with the number of `LOTS' used until all the unit cells have been used at least once on average. [Note the program Discus (Proffen & Neder, 1997) uses the same system of `LOTS'.]

Figure 5 Demonstrating the use of `LOTS' in sampling the crystal volume in calculations of diffuse scattering using the program Diffuse (Butler & Welberry, 1992). Note cyclic boundary conditions are used. (a) and (b) illustrate the use of two differently sized LOTS.

Diffuse has continued to be used in many studies since that time. The speed has increased enormously with the advent of multi-processor machines and has even been adapted to be used with very fast GPUs (graphics processing units). It has been shown (Gutmann, 2010) that the computations can be accelerated by at least one order of magnitude using modern consumer graphics cards.

4. Example of a 1D perturbed regular lattice

In this section a brief summary is given of a 1D growth disorder model that used Gaussian displacement variables to produce a 1D paracrystal (Welberry et al., 1980). This model was subsequently used in section 4.3 of Welberry (2004) to interpret the diffuse scattering in a urea inclusion compound.

In urea inclusion compounds, the urea molecules form a hydrogen-bonded network containing hexagonal channels that run along the c-direction. The channels can accommodate various kinds of long-chain molecules, of which dibromo­decane (DBD) is an interesting example. The DBD molecules pack end-to-end in an individual channel to form a pseudo-1D crystal, the repeat distance of which is not commensurate with the crystallographic repeat of the urea framework. The vertical diffuse bands in Fig. 6(a) arise from these chains of molecules. Each band is virtually uniform in intensity, because there are no correlations linking one channel with the next. It is seen that the first diffuse band is narrow, but higher-order ones get progressively broader.

Figure 6 (a) The 0kl diffraction pattern of dibromodecane/urea inclusion compound showing a sequence of progressively broader diffuse planes of scattering due to the dibromodecane molecules which form pseudo 1D crystals within the urea channels. Calculated peak profiles of 1D paracrystalline model: (b) Using only the first 10 terms in the summation over P in equation (3); (c) using all non-zero terms.

For this paracrystal model used to interpret the DBD-urea pattern the diffracted intensity consists of two parts: a Bragg intensity that comes from the average structure and a diffuse component that comes from the differences from the average.

Here, σ² is the variance of variables away from their average positions and r is the correlation between the displacements of nearest neighbouring sites. The diffuse scattering term I(k)_Diffuse is plotted in Figs. 6(b) and 6(c). It is seen that the expression for the intensity is the sum of a series of terms in increasing P. For Fig. 6(b) only the first 10 terms were used while for Fig. 6(c) all terms significantly more than zero were used. In some cases more than 100 terms were required. It should be noted that when σ² is small only the first term in the P series is significant and the shape of the diffuse peak is identical to that of an occupancy model.

This result will be referred to in §6 of the paper.

5. Modelling diffuse scattering circa 2010–2020

By 2014 when the study described in this section was published a great deal had changed from the early beginnings of our attempts to model diffuse scattering in molecular crystals described in §2. The advent of synchrotron radiation and ever improving area detectors means that fully three-dimensional diffuse scattering data can now be collected routinely for even quite small samples. Even more significant have been the advances in the computational resources that have become available. Not only are computational speeds faster but the amount of memory available has increased enormously too. The VAXstation 3100 computer mentioned in §3.4 had a maximum memory capability of 32 Mbytes whereas the laptop on which this paper is being written has 16 Gbytes of memory.

By collecting data using synchrotrons (Weber et al., 2001; Welberry et al., 2003), more reciprocal space can be mapped, both in terms of the maximum scattering vector that can be obtained and in the number of reciprocal space sections that can be obtained. In particular, use of an area detector means that whole volumes of reciprocal space can be mapped, and then any desired section obtained if required (Estermann & Steurer, 1998).

Coupled with greater computing power, this allows constructions of 3D models representing the 3D short-range ordered structure. This greatly improves modelling of highly anisotropic features. Simulations can be run for many thousands of Monte Carlo (MC) cycles, more thoroughly evaluating the efficacy of the model interactions, and the interactions themselves can be more representative, with more atom–atom interactions incorporated into the model.

Monte Carlo simulation has become the standard method of modelling, particularly for molecular crystals. This has led to the development of a standard simulation program ZMC (Goossens et al., 2011) which allows an MC model to be set up starting from a knowledge of the average structure and then used to produce an equilibrium structure with the Metropolis MC algorithm (Metropolis et al., 1953) and finally produce a calculated diffraction pattern.

5.1. Defining the model

5.1.1. Z-matrices

ZMC uses as its basis the idea that the interacting objects are molecules. Each molecule can be described using a z-matrix (an approach commonly used in ab initio molecular orbital calculations). This defines each atomic position in terms of an intramolecular coordinate set, i.e. in terms of bond lengths, bond angles and dihedral angles. Holding most of these values fixed but allowing a selected few to vary provides a mechanism for allowing molecules to show internal flexibility, while keeping other fragments of the molecule as rigid units. The molecule is then placed into its position within the model crystal by external variables which indicate which unit cell it is to go in; Cartesian coordinates x, y and z to define the position of the origin of the z-matrix; and a quaternion to define the orientation of the molecule.

5.1.2. Molecular energy

Once defined and positioned, molecules can have an internal energy dependent on the values of the internal degrees of freedom relative to their `equilibrium' positions, and an external energy dependent on the molecule's position and orientation (and therefore conformation) relative to its neighbours. Hence, the basis of ZMC is a forward Monte Carlo algorithm in which Hooke's law springs connect the molecules, and potentials are placed on selected internal degrees of freedom, usually torsional twists, within the molecule. This allows the molecules to interact with each other but also allows variations in molecular conformation to interact with molecular orientation and position. By allowing these interactions to manifest themselves through the MC algorithm a correlation structure arises in the model crystal. The final model is evaluated by calculating its diffuse scattering diffraction pattern using Diffuse (see §3.4) and comparing that with the observed scattering.

The energy of the crystal is considered to be the sum of intra- and intermolecular contributions. The intermolecular contribution to the crystal's energy comes from

where d_i is the length of vector i connecting atoms on adjacent molecules, d_0i is its equilibrium length and F_i is its force constant. The sum is over all contact vectors (cv). is the size-effect term, which allows that the equilibrium length required for the calculation may not be the average length as determined from Bragg scattering.

The intramolecular energy can be modelled using contact vectors as for intermolecular motions and by putting potentials on the internal degrees of freedom themselves.

where the G_i are the force constants for internal degrees of freedom, i and are their deviations from the equilibrium value. The G_jk are the interaction constants for interactions between the internal degrees of freedom.

5.1.3. Contact vectors

Energy calculation requires definition of the intermolecular contacts. These ought to comprise all inter-molecular (and non-bonded intra-molecular) atom–atom contacts in the system. In practice all possible atom–atom contacts up to some cut-off length may be used, but even then this would result in too many spring constants, F_i, to be used as adjustable parameters in the model. With the advent and availability of increasingly fast computing it has recently become viable to include such springs on all interatomic atom–atom vectors shorter than a conveniently chosen upper limit, e.g. 4 Å.

If many interactions are to be used in the modelling the force constants, F_i, cannot be used as the parameters to be varied in fitting the calculated diffraction pattern to the observation. Rather, adjustable parameters that define the force constants are required. One possibility is to define the spring constants using the slopes of atom–atom potentials such as the Buckingham potential

where A_jk, B_jk and C_jk are constants that depend on which types of atoms are interacting (Filippini & Gavezotti, 1993). The values for the force constants in the model, as a function of equilibrium separation, d_0i, can be determined from the derivative of equation (6). Compared to previous methods this allows the use of many more interatomic springs but far fewer adjustable parameters.

An alternative to equation (6) is the empirical formula

where R_jk is the sum of the van der Waals radii of the two atoms [taken from Bondi (1964)], d_0i is the average length of the contact vector and A, B and C are empirically determined constants. The values of A, B and C used in the simulations described here were A = 11, B = 0.4 and C = −8 since similar values have previously been refined using data from a range of molecular crystal systems (Chan et al., 2010; Chan & Goossens, 2012; Hudspeth et al., 2014).

5.1.4. Disorder in p-chloro-N-(p-chloro-benzylidene)­aniline

By combining the use of large 3D models (64 × 64 × 64 unit cells) with long MC simulations ( >> 5000 MC cycles) very high quality diffraction patterns can be obtained and the calculated intensity, I_calc, can be compared with the observed intensity, I_obs, on a pixel-by-pixel basis, excising only those pixels that an MC model of the short-range order cannot be expected to calculate (experimental artefacts and background, and the Bragg peaks themselves, for example). This allows an R factor to be calculated and minimised using the iterative automated algorithm first described by Welberry et al. (1998).

The original descriptions of the average structure of ClCl (see Scheme I) were published by Bernstein & Schmidt (1972) and Bernstein & Izak (1976). These studies showed that the molecular site in the structure was occupied by a disordered mixture of two different orientations of the molecule related by a 180° rotation about the long axis of the molecule so that the C atoms of the phenyl rings could not be clearly resolved.

Fig. 7 shows the results of automatic refinement of the structure using diffuse scattering as described in this section. It should be noted that the calculated patterns shown on the bottom two rows are of the diffuse scattering only since the average lattice was subtracted during the calculation using Diffuse. It is seen that the large crystal size coupled with the long MC simulations has produced long wavelength correlations in the structure that give TDS peaks around the Bragg positions that resemble the Bragg peaks in the observed data shown in the top row of patterns.

Figure 7 Observed and calculated diffraction patterns for two different models of p-chloro-N-(p-chloro-benzylidene)aniline. R factors are derived from a pixel-by-pixel goodness-of-fit and show that Model 2, which allows for increased molecular flexibility, is substantially better. Images reproduced from Goossens & Welberry (2014) with permission.

The results for two differing models for the interactions and molecular degrees of freedom in ClCl are shown in Fig. 7. The top row shows the observed data for three reciprocal sections, and the corresponding calculated diffuse scattering for the two Models is shown in the next two rows. Model 2 allows the molecule extra internal flexibility by decoupling the halves of the molecule, and bonding them together relatively weakly using Hooke's law contact vectors rather than covalent intramolecular bonds – effectively, treating the two halves of the molecule as separate molecules – and this is clearly crucial in improving the fit, as shown by the R factors. It is seen that viewed down the c axis the two models give almost identical agreement but viewed down a and b Model 2 is much the better.

Fig. 8(a) shows difference plots, (I_obs − I_calc), for the two models. The top half of each figure corresponds to Model 1 and the lower half to Model 2. It is interesting that for the hk0 section, although the R factors are almost equally good, the small remaining differences have quite distinct distributions implying that the data contains some remaining information that the modelling has not yet accounted for. The same is true but to an even greater extent for the other two sections (0kl) and (h0l).

Figure 8 (a) Comparison of the difference maps (Iobs − Icalc) for the two models of ClCl. Top half corresponds to Model 1 and bottom half corresponds to Model 2. (b) A scatterplot of all the atomic positions in a typical simulation superimposed on a single unit cell. Reproduced from Goossens & Welberry (2014) with permission.

Fig. 8(b) shows a scatterplot of atoms from molecules in a single one of the two molecular orientations, showing how the diffuse analysis can separate out the overlapping molecules. Here the atomic displacement parameters are all nicely shaped and even the somewhat anisotropic shape of the terminal Cl groups has been captured.

6. Calculating diffuse scattering using Scatty

The calculated diffuse scattering patterns shown in earlier sections of this paper used the program Diffuse described in §3.4. Recently Paddison wrote a paper entitled `Ultrafast calculation of diffuse scattering from atomistic models' (Paddison, 2019). This describes a new method of calculating diffuse scattering patterns from atomistic models and claims that the algorithms described can accelerate the calculations by a factor of at least 10². These claims have been substantiated by the results of three different examples described in the paper but also by a number of subsequent publications in which other authors have used the software, e.g. Roth et al. (2021) and Morgan et al. (2021). Such a large increase in speed would make the refinement of atomistic models to fit large volumes of diffuse scattering for molecular crystals much more practical. However there is one caveat: namely, that the speed of the calculation is reduced when atomic displacements are involved.

In this section we present the results of some calculations in which the performance of Scatty is compared with results previously obtained using Diffuse for an example (Thomas et al., 2007) which has been chosen because the diffuse scattering was shown to be due to some particularly large atomic displacements away from their equilibrium positions.

6.1. The Scatty program

Scatty uses the Fast Fourier Transform (FFT) algorithm (Cooley & Tukey, 1965) and also uses a fast method based on sampling theory to reduce high-frequency noise in the calculations. The program has been made freely available (Paddison, 2019).

In the original publication of Diffuse (Butler & Welberry, 1992), the use of FFT was considered and it was estimated from a simple analysis that it would be ten times faster than direct Fourier summation. However the possibility of using FFT was rejected because of the computer memory requirements which exceeded what was available in computers of the time. Now, however, typical computer memories are measured in Gbytes rather than the Mbytes that were current in 1992.

Paddison shows that the structure factor for the model crystal can be expressed in the form,

where and are Fourier transforms for which the FFT algorithm can be used:

See Paddison (2019) for full details.

This FFT-based approach used in Scatty is exact for systems in which the disorder is compositional, magnetic or in which any atomic displacements are drawn from a discrete set of values. On the other hand, if the atomic displacements can take a continuous range of values, the structure factor given by equation (8) cannot be directly evaluated by the FFT because both equations (9) and (10) contain factors of .

This problem can be addressed by expanding as a Taylor series in the same way as was done for the exponential in the 1D example of equation (3).

The FFT can only be applied to a specified order of approximation. As the order of expansion increases the number of individual FFT calculations and the amount of computer memory required increase rapidly. In the program Scatty this expansion level is set by a parameter EXPANSION_ORDER.

6.2. Diffuse scattering in pentachloronitrobenzene

In this section we present the results of some calculations in which the performance of Scatty is compared with results previously obtained using Diffuse (Thomas et al., 2007). This example has been chosen as the diffuse scattering it exhibits was shown to be due to some particularly large atomic displacements away from their equilibrium positions, as determined by Bragg scattering. This disorder has been described by Thomas et al. (2007).

Fig. 9(a) shows a plot of the average structure of pentachloronitrobenzene (C₆Cl₅NO₂, PCNB) viewed down [001] and Fig. 9(b) shows the view down [100]. The structure consists of molecular layers in which all the molecular planes lie normal to [001] and molecules in each layer are arranged on a triangular grid. The substituent groups of the molecules in the next layer along [001] are in close contact with three of the six substituent groups in the first layer. Similarly the substituent groups of the molecules in the layer below are in close contact with the other three substituent groups in the first layer. In the average structure each of the six substituent sites in every molecule contains NO ₂ and Cl.

Figure 9 The average structure of PCNB: (a) view down [001] and (b) view down [100]. Each of the six substituent sites in every molecule contains NO2 and Cl. (c) Scatterplots showing molecules in neighbouring layers along [001]. The three plots show the atom positions segregated according to whether the nearest contact involves a Cl⋯Cl, Cl⋯NO2 or NO2 contact. The graph shows the distibution of the actual distances between N⋯C, N⋯Cl and Cl⋯Cl.

Figs. 10(a) and 10(b) show the diffuse X-ray scattering patterns for the hk0 and 0kl sections, respectively. These data were used to establish a model of the disordered structure using automatic fitting of a Monte Carlo model.

Figure 10 Diffraction patterns for PCNB. (a) X-ray data for the hk0 section. (b) X-ray data for the 0kl section. (c), (d) are corresponding patterns calculated using Diffuse from a model crystal comprising (48 × 48 × 48) unit cells. (e), (f) are corresponding patterns calculated using Scatty from the same model crystal with the EXPANSION_ORDER parameter set to 9. Note in (f) c^* is vertical and b^* is horizontal.

The computer model of the disorder in PCNB that was established by Thomas et al. (2007) showed that the orientation of the molecules in each molecular site was, to a good approximation, quite random and it was concluded that the strong diffuse scattering was largely due to displacements induced by size-effect relaxations resulting from the different sizes of the substituent groups. The size of the model crystal used was 48 × 48 × 48 unit cells with three molecules per cell. This was large enough to give calculated diffraction patterns of a comparable quality to the observed patterns. The final diffraction patterns for the hk0 and 0kl sections calculated using Diffuse after 1000 MC cycles are shown in Figs. 10(c) and 10(d).

Fig. 9(c) shows scatterplots for all the atoms of two neighbouring molecules from all unit cells in the structure superposed in one unit cell. The three plots are segregated into pairs of neighbours that have close contacts between Cl⋯Cl, Cl⋯NO₂ and NO₂⋯NO₂, respectively. It is seen that, while for the Cl–Cl contact the planes of the two molecules lie fairly parallel to the (001) plane, for the NO₂⋯NO₂ contact the planes are tilted to significantly increase the average NO₂⋯NO₂ distance. The graph in Fig. 9(c) shows the distribution of actual distances between N⋯N, N⋯Cl and Cl⋯Cl. The distances involved in the interaction are summarised in Table 1.

Table 1
The effect of size on atom–atom vectors comprising the intermolecular contact shown in Fig. 10(c)

Vector type	Average structure (Å)	Mean (Å)	σ	Δd	% change
N⋯N	3.817	4.404	0.27	0.587	+15
N⋯Cl	3.919	4.201	0.27	0.282	+7
Cl⋯N	3.919	4.201	0.27	0.282	+7
Cl⋯Cl	4.035	3.951	0.27	−0.084	−2

6.3. Diffraction patterns of PCNB computed with Scatty

The final atomic coordinates output from the MC simulation of PCNB were used as input to the program Scatty. This was quite straightforward. The calculations were made on a MacBook Pro computer with a 2.9 GHz Intel Core i7 processor and 16 GBytes of memory.

Because of the large atomic displacements that were present in the PCNB model, diffraction patterns for the hk0 and 0kl sections were computed with a range of different values of the EXPANSION_ORDER parameter. Since the computation time increases with the order number, first examples used a small expansion number and this was gradually increased. For the hk0 section it was found that it was necessary to go to a value of 9 before the resulting pattern became at all comparable to the pattern obtained with Diffuse. The Scatty pattern is shown in Fig. 10(e). Even so, the pattern rapidly deteriorates outside the Q-range shown, an indication that level 9 is the bare minimum that is required. At this level of calculation the number of products in the Taylor expansion of displacements was 219, the amount of memory required was 45 Gbytes and the total time for the calculation was 18.5 min. (see Table 2). This compares with a time of 34 min for the comparable Diffuse calculation [Fig. 10(c)].

Table 2
Summary of the different Scatty computations for the examples shown in Figs. 10 and 11

The columns show: the model crystal size; the order of the exponential expansion; the amount of memory used; the number of products involved in the exponential expansion; the time used in calculating the FFTs and the total CPU usage in seconds.

Fig.	Size	Order	GByte	No. Prod.	FFT	Total CPU
10(e)	48 × 48 × 48	9	45.5	219	623	1109
10(f)	48 × 48 × 48	9	45.4	219	621	5190
11(a)	24 × 24 × 24	18	34.5	1329	277	6736
11(b)	24 × 24 × 24	9	5.7	219	35	417

The situation turned out to be not as good for the hk0 section. This appears to be due to the fact that it is in this section that the large size-effect displacements feature most prominently. Fig. 10(f) shows the computed pattern for this section using level 9 expansion. Here it is seen that close to the horizontal ( b^*) axis the pattern looks quite good but along c^* the error in the calculation pattern becomes increasingly bad. Beyond a certain point the errors become too large and Scatty produces no output (grey shading). The computation time for this section was 86 min (see Table 2).

In fact it was found that with the current memory available in the MacBook Pro computer, a value of EXPANSION_ORDER = 9 was the maximum that could be achieved using the (48 × 48 × 48) unit-cell crystal. As an experiment therefore some calculations were carried out with a much smaller crystal (a 24 × 24 × 24 section of the larger crystal was used). It was then found possible to use EXPANSION_ORDER = 18 for 0kl. The result of this calculation is shown in Fig. 11(a). Even with this large value of the expansion order the valid calculation does not extend to fill the circle completely. However the diffuse blob around the (0012) position is revealed and it is seen this is diminished in intensity relative to that around (009) in agreement with the observed pattern.

Figure 11 Attempts to improve the calculation of the 0kl section. (a) Using a smaller crystal (24 × 24 × 24) with EXPANSION_ORDER = 18. (b) Using a (24 × 24 × 24) crystal with EXPANSION_ORDER = 9 but with all c^* components of the atom displacements (away from their average positions) reduced by a factor of 1.5. Note c^* is vertical and b^* is horizontal.

A second experiment was carried out in which all of the atomic displacements in the c^* direction away from their average sites were reduced by a factor of 1.5. A calculation of the 0kl section for this modified crystal using level 9 expansion is shown in Fig. 11(b). Now, with the reduced displacement magnitudes the calculated pattern extends to fill the whole circle. Moreover the amount of memory required is much reduced and the computation time was only 7 min. What is noticeable, though, is that the magnitude of the intensities of the large diffuse blobs around the (003), (006), (009), (0012), (0015) positions do not decay towards higher Q. It is seen that in the X-ray pattern and the Diffuse calculation that the (0012) blob is very weak. This drop in intensity at high Q is indicative of the magnitude of the size-effect distortion and was crucial in the original disorder model refinement.

7. Conclusion

In this paper a brief account has been given of how our modelling of diffuse scattering from molecular crystals has developed over the last 40+ years. Not only has the advent of synchrotron radiation and modern area detectors enabled the measurement of high quality three-dimensional diffuse scattering data but the accompanying exponential growth of computational power, that has become available for processing and analysing it, has had an enormous impact.

The example described in §5.1.4 shows that the goal of trying to make a full 3D refinement of a model for the disorder in a molecular structure is close to being realised. As the simulation methodology improves it might be expected that finer and finer detail may become accessible, just as occurred in conventional crystallography some 60–70 years ago. The analysis of diffuse scattering thus appears to be poised to enter a new period of expansion and exploitation, analogous to that which occurred for Bragg scattering 60 to 70 years ago.

The exciting new method of calculating diffuse scattering using the program Scatty will undoubtedly have a major role to play into the future. The issues that have been raised by the example discussed in §6.3 regarding the performance of Scatty when large atomic displacements are present do not detract from the overall promise of the new method. For many systems the problems will not apply but for others it is likely that with the possibility of faster processors, the use of parallel processing and larger amounts of available memory the issues will be less significant. However it is important that current users be aware of the issue.