Ultrafast Calculation of Diffuse Scattering from Atomistic Models

# Ultrafast Calculation of Diffuse Scattering from Atomistic Models

## Abstract

Diffuse scattering is a rich source of information about disorder in crystalline materials, which can be modelled using atomistic techniques such as Monte Carlo and molecular dynamics simulations. Modern X-ray and neutron scattering instruments can rapidly measure large volumes of diffuse-scattering data. Unfortunately, current algorithms for atomistic diffuse-scattering calculations are too slow to model large data sets completely, because the fast Fourier transform (FFT) algorithm has long been considered unsuitable for such calculations Butler_1992 (). Here, I present a new approach for ultrafast calculation of atomistic diffuse-scattering patterns. I show that the FFT can be actually be used to perform such calculations rapidly, and that a fast method based on sampling theory can be used to reduce high-frequency noise in the calculations. I benchmark these algorithms using realistic examples of compositional, magnetic, and displacive disorder. They accelerate the calculations by a factor of at least , making refinement of atomistic models to large diffuse-scattering volumes practical.

## I Introduction

Disorder plays an increasingly important role in our understanding of crystalline materials. Whereas conventional crystallography is primarily concerned with the average positions of atoms or molecules, local deviations from the average structure are fundamental to the properties of many important systems. Topical examples include fast-ion conductions Keen_2002 (), frustrated magnets Fennell_2009 (); Paddison_2013a (), the polar nanodomains of lead-based perovskite ferroelectrics Xu_2004 (); Pasciak_2012 (), and the orbital correlations of colossal magnetoresistance manganites Shimomura_1999 (). Diffuse scattering—the weak features observed beneath and between the Bragg peaks in scattering experiments—plays a central role in understanding structural disorder Billinge_2007 (). Whereas Bragg peaks arise from the ideal periodicity of the average structure and yield information about single-particle correlations, diffuse scattering arises from local deviations from the average structure and yields information about pairwise correlations. In disordered crystals, strong pairwise correlations are often present at the nanoscale, yielding highly-structured diffuse-scattering patterns Keen_2015 (); Overy_2016 (); Weber_2012 (). Disorder in crystalline materials can often be divided into three broad categories: local atomic displacements, which may be either static or dynamic (phononic); local variations in chemical composition or site occupancy such as site mixing and atomic vacancies; and disordered spin arrangements in magnetic materials. In all cases, the modulation of the diffuse-scattering intensity provides vital information about how the relevant degrees of freedom—whether displacive, compositional, or magnetic—are locally correlated.

Recent advances in instrumentation have allowed rapid collection of large three-dimensional (3D) diffuse-scattering datasets from single crystals. In neutron-scattering experiments, such datasets can be measured in a matter of days using time-of-flight instruments with large detector coverage, such as SXD at the ISIS Neutron Source Keen_2006 () and Corelli at ORNL Ye_2018 (). In X-ray scattering experiments, area detectors such as the Pilatus Kraft_2009 () allow such datasets to be obtained in minutes. These experimental developments have been coupled with advances in 3D atomistic modelling. Atomistic models may be generated using Monte Carlo or molecular dynamics simulations based on a set of interaction parameters, which are determined either from first-principles simulations Gutmann_2013 () or by fitting to experimental data Welberry_1998 (); Welberry_2008 (); Weber_2002 (). Alternatively, atomistic models may be generated by fitting atomic positions or magnetic-moment orientations directly to experimental data, in an approach called reverse Monte Carlo refinement McGreevy_1988 (); Tucker_2007 (); Paddison_2012 (). Improvements in computing power have greatly reduced the time required to generate atomistic models using these approaches.

A key step in atomistic diffuse-scattering analysis is calculating the diffuse-scattering intensity, which must usually be repeated several times as the model is improved. In principle, this involves the straightforward procedure of taking the Fourier transform of a set of atomic positions. In practice, however, such calculations are limitingly slow, because of two main problems. The first problem is that the fast Fourier transform (FFT) algorithm—which can yield enormous increases in speed compared to direct Fourier summation—has been considered unsuitable for these calculations Butler_1992 (). This is because the FFT requires sampling atomic positions on an equally-spaced grid Cooley_1965 (), but a grid that sampled local atomic displacements accurately would need so many samples that the computational advantages of the FFT would be negated Butler_1992 (); Neder_1989 (); Rahman_1989 (). The second problem is that diffuse-scattering calculations are often marred by high-frequency noise. This occurs because the entire atomistic model is assumed to scatter coherently, which often allows the many essentially-uncorrelated atomic displacements at large interatomic distances to dominate the calculation Butler_1992 (). Current software addresses this problem by dividing the supercell into many smaller regions called “sub-boxes” or “lots” and averaging the scattering intensity over sub-boxes Butler_1992 (); Proffen_1997 (). This approach is effective at reducing noise, but has a large computational cost that can only be reduced by parallelization Gutmann_2010 (). Because of their computational cost, diffuse-scattering calculations have typically been restricted to a few 2D sections of reciprocal space Welberry_2008 (). In most cases, therefore, the information content of 3D diffuse-scattering data has yet to be fully explored.

Here, I present an approach to diffuse-scattering calculations that seeks to address the computational limitations of current algorithms. I present two main results. First, I critically reassess the assumption that the FFT is not useful for atomistic diffuse-scattering calculations. I show that, on the contrary, the FFT can be used to perform such calculations rapidly. This FFT approach is exact for systems in which the disorder is compositional, magnetic, or in which atomic displacements are drawn from a discrete set of values; if atomic displacements can take a continuous range of values, then the FFT can be applied to a specified order of approximation. Second, I show that the desirable noise-reduction properties of the “sub-box” approach can also be obtained by a faster method based on sampling theory. I demonstrate the practicality of these algorithms by presenting example calculations based on real materials, and by providing a program Scatty for ultrafast diffuse-scattering calculations (available from the author on request). These developments accelerate atomistic diffuse-scattering calculations by several orders of magnitude compared to current algorithms, making refinement of atomistic models to 3D diffuse-scattering data a practical possibility.

This article is structured as follows. In Section II, I review the scattering equations and aspects of sampling theory in the context of atomistic models. In Section III, I show how the FFT can be used to accelerate atomistic diffuse-scattering calculations. In Section IV, I present example calculations for model systems that exhibit occupational, magnetic, and displacive disorder, in which highly-structured diffuse scattering is driven by “ice rules”. In Section V, I show how noise in the calculations can be reduced using a resampling approach. I conclude in Section VI with a discussion of the implications of this work.

## Ii Atomistic simulations

In atomistic simulations, disorder is modelled using a “virtual crystal”. This is a supercell of the crystallographic unit cell, typically containing atomic positions, in which each atomic position is decorated by variables corresponding to its occupancy by a certain element, its displacement from its average position, and/or the orientation of its magnetic moment. I consider a supercell that consists of crystallographic unit cells, where , , and are the numbers of unit cells parallel to the , and crystal axes, respectively. An atomic position in the supercell is given by

 r=n1r1a+n2r2b+n3r3c;0≤rα<1. (1)

An interatomic vector is denoted , with components . A general wavevector is given by

 Q =ha∗+kb∗+lc∗ (2) ≡Q1n1a∗+Q2n2b∗+Q3n3c∗;Qα∈R, (3)

where reciprocal-lattice vectors are defined as , etc., where is the unit-cell volume. Atomistic simulations usually impose periodic boundary conditions to avoid edge effects, and I assume throughout this article that this is the case. The Bragg positions of the periodic supercell are then given by

 G =G1n1a∗+G2n2b∗+G3n3c∗;Gα∈Z, (4)

and become increasingly closely spaced as the number of unit cells in the supercell is increased.

The energy-integrated coherent neutron-scattering intensity from a real crystal is given, in the kinematic approximation, by

 I(Q) =1N∑j,k⟨bjbkexp[iQ⋅(rj−rk)]⟩ (5) ≡∑Δrρ(Δr)exp(iQ⋅Δr), (6)

where is the coherent neutron scattering length of atom , and is the 3D pair-distribution function Weber_2012 (). Angle brackets denote spatial and temporal averaging, which can be approximated in atomistic modelling by averaging over many supercells. For X-ray scattering, the are replaced by -dependent X-ray form factors, . By factorising the double summation in equation (5), the scattering intensity can be rewritten as

 I(Q)=1N⟨|F(Q)|2⟩, (7)

where the structure factor is a discrete Fourier transform,

 F(Q)=N∑j=1bjexp(iQ⋅rj). (8)

In general, the scattering intensity from a disordered crystal can be separated into its Bragg and diffuse contributions,

 I(Q) =IBragg(Q)+Idiffuse(Q), (9) =1N|⟨F(Q)⟩|2+1N⟨|F(Q)−⟨F(Q)⟩|2⟩, (10)

which arise from the average structure and local modulations away from the average, respectively Frey_1997 ().

Atomistic models inevitably contain far fewer unit cells () than do real crystals (), and it is important to consider how this affects scattering calculations. Fundamentally, a periodic supercell only contains information about pair correlations for which , where . Scattering calculations may account for this in two different ways. First, if the scattering is calculated using equation (5), then the separation between pairs is taken to be the shortest separation between their periodic images; this is called the nearest-image convention and involves replacing with where necessary to ensure . This procedure can be conceptualized as multiplying the pair-distribution function for an infinite tiling of the supercell, , by a cutoff function, so that

 ρ(Δrα)=ρinf(Δrα)⋅cutoff(Δrα), (11)

where is a rectangular function equal to for and zero elsewhere. In this case, the calculated is a smooth function of wavevector, which contains finite-size artifacts because the contribution of the average structure to is truncated. In contrast, is free from finite-size artifacts provided that the local modulations are short-ranged compared to . The second approach is to calculate the scattering using equations (7) and (8). In this case, the nearest-image convention cannot be applied. Consequently, the scattering intensity may only be sampled at the supercell Bragg positions , for which is unchanged by the nearest-image convention because . The two approaches become equivalent for real crystals, which are large enough that the boundary conditions become irrelevant.

In practice, it is much faster to calculate the scattering intensity using equations (7) and (8)—which require only a single summation over atomic positions—than to use equation (5), which requires a double summation. However, the resulting limitation that only may be calculated is often too restrictive, especially when a direct comparison with experimental data is required Butler_1992 (). To solve this problem, I consider the information content of the scattering pattern. Applying the convolution theorem to equation (11) yields

 I(Q) =∞∑G=−∞I(G)Wsinc(Q−G), (12)

where the weight function

 Wsinc(Q−G) =8r3cut3∏α=1sinc[2πrcut(Qα−Gα)] (13)

is the Fourier transform of the rectangular function, is the Fourier transform of sampled at , and Whittaker_1915 (); Shannon_1949 (). In the context of sampling theory, equations (12) and (13) are known as the Whittaker-Shannon interpolation formulae. They show that the scattering intensity can, in principle, be reconstructed at any wavevector , given only its samples at the supercell Bragg positions . This is possible because the supercell Bragg positions sample at its Nyquist rate Sayre_1952 ().

In the rest of this article, I will first show how can be calculated rapidly using the FFT algorithm. I will then show how can be estimated in practice, by modifying equations (12) and (13).

## Iii Fast Fourier transform

### iii.1 Requirements for the FFT

The FFT is an algorithm that calculates the discrete Fourier transform rapidly Cooley_1965 (); Singleton_1969 (). In the 3D case relevant for crystalline materials, it calculates a function of the form

 Xk =∑RxRexp(ik⋅R). (14)

Unlike the discrete Fourier transform, for which the position vector and wavevector can take any real values, the FFT imposes two restrictions. First, the must be arranged on an equally-spaced grid, so that

 R=R1a+R2b+R3c; Rα∈{0,1,...,Nα−1}, (15)

where denotes the number of grid points parallel to crystal axis in the supercell. Second, the values of at which the Fourier coefficients may be evaluated are given by

 k=k1N1a∗+k2N2b∗+k3N3c∗;kα∈{0,1,...,Nα−1}. (16)

Below, I will identify the grid points with lattice points, so that and each grid point is the origin of a particular unit cell within the supercell. In this case, the values of given by equation (16) are a subset of the supercell Bragg positions that lie within a fundamental domain of the reciprocal lattice of the supercell.

The computational cost of the discrete Fourier transform is on the order of , whereas for the FFT it is on the order of —a huge saving for large . Unfortunately, the atomic positions in a crystal do not generally satisfy equation (15), except in the special case where atoms occupy a Bravais lattice without disorder. In early work, it was assumed that this restriction could only be addressed by using a fine grid encompassing all the atomic positions to some specified accuracy; however, the number of grid points required in 3D was found to be be prohibitively large Butler_1992 (); Neder_1989 (); Rahman_1989 (). Consequently, the FFT has been considered unsuitable for atomistic diffuse-scattering calculations, and is not implemented in current software such as Diffuse Butler_1992 (), Discus Proffen_1997 (), ZMC Goossens_2011 (), and Zods Michels-Clark_2014 ().

I now show that the scattering intensity can actually be rewritten in a form that allows the FFT to be efficiently applied, by making use of the underlying periodicity of the average structure. To emphasise this periodicity, I write the position of an atom in the supercell as

 rR,μ,i≡R+rμ+uR,μ,i, (17)

where henceforth denotes a lattice point given by equation (15) with ; is the average position of site within the crystallographic unit cell; and is the local displacement of the atom of element belonging to site at lattice point . A Bragg position of the supercell can be written

 G≡H+k, (18)

where henceforth denotes a wavevector given by equation (16) with , and

 H=H1a∗+H2b∗+H3c∗;hα∈Z (19)

is a Bragg position of the crystallographic unit cell. By substituting equations (17) and (18) into equation (8), and using the fact that , the structure factor can be expressed as

 F(G)=∑R,μ,iδR,μ,ibμ,iexp[iG⋅(rμ+uR,μ,i)]exp(ik⋅R), (20)

where is equal to if site at lattice point is occupied by an atom of element , and is otherwise zero, and is a coherent neutron-scattering length. I also define the difference between the local occupancy and the average occupancy,

 aR,μ,i =δR,μ,i−cμ,icμ,i, (21)

where is the average occupancy of site by atoms of element . This separation of the structure into an average part and a local modulation will allow the Bragg and diffuse contributions to the strucure factor to be separated. From equations (20) and (21), the structure factor is given by

 F(G)=∑μ,i[Uk,μ,i(G)+Ak,μ,i(G)]cμ,ibμ,iexp(iG⋅rμ), (22)

where I have defined a pair of Fourier transforms for each site and element ,

 Uk,μ,i(G) =∑Rexp(iG⋅uR,μ,i)exp(ik⋅R), (23) Ak,μ,i(G) =∑RaR,μ,iexp(iG⋅uR,μ,i)exp(ik⋅R), (24)

which are also used in the “modulation wave” approach to diffuse-scattering analysis Withers_2015 (). The Bragg structure factor is given by

 ⟨F(G)⟩=n1n2n3δG,H∑μ,iTμ,i(G)cμ,ibμ,iexp(iG⋅rμ), (25)

where

 Tμ,i(G)=1n1n2n3∑R⟨exp(iG⋅uR,μ,i)⟩ (26)

is the Debye-Waller factor. I will show below how the FFT can be applied to calculate equation (22).

### iii.2 Compositional/occupational disorder

I consider first a scenario in which the disorder is purely compositional or occupational (these terms will be used interchangeably), so that all atomic displacements can be set to zero. In this case, the structure factor given by equation (22) reduces to

 Fc(G) =⟨Fc(G)⟩+∑μ,iak,μ,icμ,ibμ,iexp(iG⋅rμ), (27)

in which

 ak,μ,i=∑RaR,μ,iexp(ik⋅R) (28)

is given by equation (24) with , and is given by equation (25) with . The diffuse intensity from compositional disorder is then given by

 Ic,diffuse(G)=1N⟨|Fc(G)−⟨Fc(G)⟩|2⟩, (29)

where angle brackets denote averaging over supercells. More effective averaging is obtained by using many smaller supercells instead of a few larger ones; I recall that the former approach does not produce finite-size artifacts in diffuse-scattering calculations, provided that for the supercell exceeds the correlation length of the local modulations.

A key observation of this article is that equation (28) has the same form as equation (14), with the substitutions and . Consequently, equation (28) can be readily calculated by the FFT. This is possible because each site in the unit cell forms a Bravais lattice; hence, equation (28) is periodic in reciprocal space Borie_1971 (); Hayakawa_1975 (). Equation (28) can also be applied to models in which local atomic displacements are present, provided they can take only relatively few different magnitudes and directions. In such cases, the number of atoms in the crystallographic unit cell is increased according to the number of possible displacements, and each displaced atomic position is assigned to a site in the crystallographic unit cell; in this way, displacive disorder is mapped to occupational disorder in a “split-site” model. Because many atomistic models of real materials assume discrete displacement distributions Welberry_2008 (), equation (28) is often applicable in practice.

The computer time required to calculate at wavevectors using the FFT is approximately proportional to the sum of two terms, where and are the times required to calculate equations (27) and (28), respectively, and is the number of atoms in the unit cell. Application of the discrete Fourier transform to equation (20) requires a time . In typical simulations, , , and may range from for scattering planes to more than for large volumes. In practice, I will show in section IV that the FFT accelerates typical scattering calculations by two to three orders of magnitude compared to the discrete Fourier transform.

### iii.3 Compositional/occupational and displacive disorder

If atomic displacements can take many possible magnitudes or directions, the structure factor given by equation (22) cannot be directly evaluated by the FFT, because and contain factors of and hence do not have the same form as equation (14). I now show that this problem can be addressed by expanding as a Taylor series,

 exp(iG⋅uR,μ,i)=∞∑n=0inn!(G⋅uR,μ,i)n. (30)

Writing the scalar product in terms of its components yields

 G⋅uR,μ,i =3∑α=1GαuαR,μ,i, (31) (G⋅uR,μ,i)2 =3∑α1=13∑α2=1Gα1Gα2uα1R,μ,iuα2R,μ,i, (32) (G⋅uR,μ,i)3 =3∑α1=13∑α2=13∑α3=1Gα1Gα2Gα3uα1R,μ,iuα2R,μ,iuα3R,μ,i, (33)

and so on; evidently, the term of order in the Taylor expansion contains a sum of all products of components of and components of . Crucially, the products depend only on atomic position, and their Fourier transforms can therefore be evaluated by the FFT. The generalisation of equations (31)–(33) for the term of order is

 (G⋅uR,μ,i)n =∑α(n)Gα(n)Uα(n)R,μ,i, (34)

where , and , and . Substituting equations (30) and (34) into equations (23) and (24), and the results into equation (22), I obtain the structure factor

 (35)

where I have defined the Fourier transforms

 Uα(n)k,μ,i =∑RUα(n)R,μ,iexp(ik⋅R), (36) Aα(n)k,μ,i =∑RaR,μ,iUα(n)R,μ,iexp(ik⋅R). (37)

Equations (36) and (37) have the same form as equation (14) and can therefore be evaluated by the FFT. The Bragg structure factor is given by

 ⟨F(G)⟩ (38)

where only even- terms are included in the Taylor expansion because the odd- terms average to zero. Finally, the diffuse intensity is obtained as

 Idiffuse(G)=1N⟨|F(G)−⟨F(G)⟩|2⟩. (39)

In practice, it is necessary to truncate the Taylor expansion of at a finite order of approximation. The number of FFTs that must be performed to evaluate the term of order is equal to , which increases rapidly with . Nevertheless, because of the favourable scaling of the FFT, calculations with negligible loss of accuracy still allow large performance improvements over the discrete Fourier transform, as I will show in section IV.1

### iii.4 Magnetic disorder

I now consider a magnetic system in which magnetic moments decorate atomic positions. Because the neutron has a magnetic moment, neutron scattering is a powerful technique to study local magnetic correlations in materials such as spin liquids, spin glasses, and frustrated magnets. The magnetic structure factor for neutron scattering is a vector quantity Squires_1978 (); Lovesey_1987 (),

 Fmag(G)=∑μ,ifmagμ,i(|G|)Tμ,i(G)Mk,μ,iexp(iG⋅rμ), (40)

where is the magnetic form factor for atoms of element belonging to site , and coupling between magnetic and displacive variables is included only via the Debye-Waller factor . The Fourier transform of the magnetic moments is given by

 Mk,μ,i=∑RMR,μ,iexp(ik⋅R). (41)

It is straightforward to calculate equation (41) using the FFT, as has been noted previously Conlon_2010 (); the only difference in implementation between compositional and magnetic cases is that the FFT is applied for each component of the vector in the magnetic case. The magnetic neutron-scattering intensity is given by

 Imag(G)=CN⟨∣∣F⊥mag(G)∣∣2⟩, (42)

where the constant  barn, and

 F⊥mag=^G×Fmag×^G (43)

is the projection of the magnetic structure factor perpendicular to . If the magnetic moments do not show long-range order, then the magnetic scattering will be entirely diffuse; this is the case in, e.g., spin-ice materials Bramwell_2001 (); Fennell_2009 ().

## Iv Example calculations

To benchmark the FFT-based approach, I present diffuse-scattering calculations for simple examples of occupational, magnetic, and displacive disorder. In all cases, I use the same underlying model of “ice rules” on a pyrochlore network of corner-sharing tetrahedra.

As an example of occupational disorder, I consider an idealized model of proton disorder in cubic water ice Pauling_1935 (). Oxygen atoms are ordered and lie at the centres of the tetrahedra. Protons are disordered according to the “ice rule” that each oxygen has two covalently-bonded protons close to it, such that these two O—H bonds point towards vertices of the tetrahedron. The model is disordered because there are six ways of satisfying this constraint on each tetrahedron, which gives rise to the well-known zero-point entropy of water ice Pauling_1935 (); Giauque_1936 (). The disorder is occupational because tetrahedra have four possible proton sites, of which only two are occupied for any given tetrahedron.

As an example of magnetic disorder, I consider the spin-ice state observed in materials such as HoTiO and DyTiO, in which magnetic Ho or Dy ions occupy a pyrochlore network. The Ho/Dy magnetic moments are constrained by the crystalline electric field to point towards or away from the tetrahedra centres, and the magnetic interactions are such that, at low temperatures, two moments on each tetrahedron point towards its centre and two point away Bramwell_2001 (). The resulting degeneracy of spin arrangements leads to a zero-point entropy equivalent to that of water ice Ramirez_1999 ().

Finally, as an example of displacive disorder, I consider the pyrochlore material YMoO Gardner_1999 (), in which an “orbital ice” state was recently proposed. In the proposed model, Jahn-Teller-active Mo ions locally displace parallel to Mo—Mo vectors such that two Mo ions on each tetrahedron displace towards each other Thygesen_2017 (). The resulting degeneracy is again equivalent to water ice and spin ice. The ice-rules models are all characterized by the selection of a single edge on each tetrahedron, which may be identified with the two protons close to the oxygen in water ice, the two spins pointing away from the tetrahedron centre in water ice, or the two Mo ions that displace towards each other in YMoO Keen_2015 ().

The diffuse scattering from each model was simulated as follows. First, cubic supercells that obeyed the ice rule were generated using a Metropolis Monte Carlo algorithm. Each supercell contained conventional unit cells ( vertices of the pyrochlore network). Second, these supercells were mapped to the disorder models discussed above. For the water-ice model, I used a lattice constant  Å, a covalent O—H bond length of  Å, and the neutron-scattering length for H. For the spin-ice model, I used  Å and the neutron magnetic form factor for Ho. For the orbital-ice model, I used  Å, a Mo displacement magnitude of  Å Thygesen_2017 (), and the atomic X-ray form factor for Mo. Finally, the diffuse intensity was calculated in the scattering plane on a grid of pixels over the reciprocal-space range and , which contained supercell Bragg positions. The diffuse intensity was calculated for the 100 supercells separately, and these results averaged to improve the statistical quality of the final calculation. The FFT was implemented using an updated version of the Fortran code by Singleton Singleton_1969 (). For FFT calculations, equations (27)–(29) were used for the water-ice model, and equations (40)–(43) for the spin-ice model. For the orbital-ice model, the diffuse intensity was calculated in two different ways using the FFT: first, by mapping the displacive disorder to occupational disorder and applying the exact equations (27)–(29); and second, by applying a fifth-order Taylor approximation using equations (35)–(39). For comparison, the discrete Fourier transform was also used to calculate equations (20) and (40) without taking advantage of reciprocal-space periodicity. All calculations were performed on a laptop with a  GHz Intel Core i5 processor, and the CPU time measured.

Figure 1(a)i–(c)i shows water-ice, spin-ice, and orbital-ice scattering calculations, respectively. The results for water ice and spin ice are in good agreement with published calculations Wehinger_2014 (); Fennell_2009 (). Calculations using the discrete Fourier transform and the FFT gave identical results, as required. Crucially, however, the FFT calculations are much faster; e.g., the CPU time required for the water-ice calculation was  s for the FFT compared to  s for the discrete Fourier transform—an improvement of two orders of magnitude. Moreover, although a total of 55 FFTs were required for the fifth-order approximation in the orbital-ice calculation, a reduction in CPU time of two orders of magnitude was still obtained. The largest error in any pixel for this fifth-order approximation was less than , confirming that the Taylor expansion can be used for quantitative work Butler_1993 (). I emphasize that, while Figure 1 shows a scattering plane, large scattering volumes can be calculated in a few minutes; e.g., for the water-ice model, a cube in reciprocal space extending from and containing supercell Bragg positions was calculated in  s.

## V Noise reduction

### v.1 Method of sub-boxes

The scattering patterns shown in figure 1(a)i–(c)i, while of good quality, are affected by two problems. First, because the diffuse scattering has been calculated only at supercell Bragg positions, it is not a smooth function of wavevector. Second, the calculations visibly affected by high-frequency noise (“speckle”), even though the results shown have been averaged over 100 independent scattering calculations to reduce noise. As discussed in section I, this noise occurs because the calculation can be dominated by the many essentially-uncorrelated atom pairs that are large distances apart.

A proven approach to address these problems is to divide the supercell into a set of smaller overlapping regions, called “sub-boxes” or “lots” Butler_1992 (); Proffen_1997 (). Each sub-box has open boundary conditions, and contains unit cells parallel to each crystal axis ( unit cells in total). The value of should exceed the correlation length of the local modulations that generate the diffuse scattering, but be smaller than or equal to , in order to avoid contributions from periodic images of the supercell. The intensity for each sub-box is calculated separately using the discrete Fourier transform, and the results averaged to estimate . In the original algorithm, sub-boxes are distributed at random within the supercell, and the number of sub-boxes is chosen so that each atom is included at least once Butler_1992 (). However, this random-sampling approach has the disadvantage that a different answer is obtained each time the calculation is run. A recent algorithm addressed this problem by using every lattice point in the supercell as the origin of a sub-box, while avoiding unnecessary repetition of phase-factor calculations Paddison_2013b (). I used this algorithm to perform scattering calculations for the example models discussed in section IV, using a sub-box length unit cells.

Example scattering calculations using the method of sub-boxes are shown in figure 1(a)ii–(c)ii. The use of open boundary conditions for each sub-box means that the calculation is not restricted to the Bragg positions of the supercell. Moreover, the calculation is essentially free of visible noise. There are two reasons for these successes. First, it is possible in general to choose , which excludes longer-range correlations from the calculation Butler_1992 (). However, since I take , this reason does not apply in the present example. In fact, the method of sub-boxes is effective here because each sub-box has open boundary conditions. To understand this, I recall that for a periodic supercell, the number of distinct pairs of lattice points separated by unit cells (parallel to crystal axis ) is equal to ; the cutoff function in equation (11) is therefore a rectangular function, as shown in figure 1(d)i. In contrast, if the supercell has open boundary conditions, then the number of distinct pairs of lattice points separated by unit cells is equal to . Consequently, the cutoff function for an open supercell is a triangular function that is equal to unity at and zero at . Hence, in the method of sub-boxes, the cutoff function is also triangular and decays to zero at , as shown in figure 1(d)ii. This is the main result of this section. It explains how the method of sub-boxes reduces noise, because as the separation of atom pairs increases, their pair correlations are suppressed in the scattering calculation.

The analysis above also reveals two limitations of the method of sub-boxes. First, it is very much slower than the FFT. Indeed, it is also much slower than the discrete Fourier transform used to calculate (figure 1(a)i–(c)i): this is mainly because the grid of wavevectors samples reciprocal space at finer intervals than the Nyquist rate defined by the spacing of the supercell Bragg positions, and hence calculates redundant information. Second, its effectiveness at reducing noise has the trade-off that the scattering is artificially blurred. Although this blurring can be reduced by increasing , it is never entirely absent, because the triangular cutoff function is never flat. This limitation may be significant if accurate estimates of the correlation length are required.

### v.2 Lanczos resampling

I now develop an approach that addresses the limitations of the method of sub-boxes, while maintaining its noise-reduction properties. Whereas the method of sub-boxes works directly in real space, I will show that its effects can be replicated efficiently by applying a resampling filter in reciprocal space.

I first recall that can, in principle, be calculated at any wavevector given only its samples at the supercell Bragg positions , by applying the Whittaker-Shannon interpolation formulae, equations (12) and (13). Consequently, can be calculated using the FFT and then resampled to obtain on a fine grid. An apparent problem with this resampling approach is that equation (12) requires summation over an infinite number of supercell Bragg positions to reproduce the discontinuity in the rectangular cutoff function at . If the infinite summation is replaced by a finite summation, then large truncation artifacts (Fourier ripples) are introduced in the cutoff function, which do not vanish even as the number of summed terms becomes very large; this undesirable effect is known as the Gibbs phenomenon Gibbs_1899 (). However, the key result of section V.1—that the method of sub-boxes in effect applies a triangular cutoff function—suggests that a rectangular cutoff function is not necessary or even desirable for diffuse-scattering calculations. The triangular cutoff function, and other suitably-chosen continuous functions, can be approximated by using a finite summation in equation (12). Consequently, the scattering intensity can be approximated as

 I(Q) ≈∑GI(G)W(Q−G)∑GW(Q−G), (44)

where the weight function is the Fourier transform of the cutoff function, and the summation runs over supercell Bragg positions in the vicinity of . This is the basis for the technique called Lanczos resampling Lanczos_1956 (); Duchon_1979 (), which I discuss below.

The effects of Lanczos resampling in real and reciprocal space can be derived in three steps Duchon_1979 (). First, the desired form of the cutoff function is specified. Second, is determined as the Fourier transform of this idealized cutoff function. Third, is back-Fourier-transformed numerically to determine the actual cutoff function obtained for a finite summation range. Following Lanczos Lanczos_1956 (), I consider an idealized cutoff function that is equal to unity for and decays linearly to zero over the range , where is an integer Duchon_1979 (). Hereafter, the nominal cutoff separation is the separation at which the idealized cutoff function equals . Since this cutoff function is the convolution of two rectangular functions, its Fourier transform is the product of two sinc functions, so that

 W(Q−G)=Wsinc(Q−G)L(Q−G), (45)

where is given by equation (13) with , and

 L(Q−G) ={∏αsinc[π(Qα−Gα)/m]if |Qα−Gα|

is called the window function. Because equation (46) is set to zero outside the central lobe of its sinc function Duchon_1979 (), the summation in equation (44) is only taken over the values of for which , where square brackets denote the floor function. This approach is called windowed-sinc filtering. It has been widely used for image filtering and digital-signal processing Getreuer_2011 (); Marschner_1994 (). It was applied to diffuse-scattering calculations for the first time by the present author Paddison_2018 (); a more detailed analysis is given here.

Example scattering calculations using Lanczos resampling with are shown in figure 1(a)iii–(c)iii and calculations with are shown in figure 1(a)iv–(c)iv. In each case, was calculated using the FFT, and was then estimated using equations (13) and (44)–(46). The corresponding cutoff functions are shown in figure 1(d)iii and 1(d)iv. The idealized cutoff function for the case is a triangular function, and the actual cutoff function resembles a smoothed version of it. Consequently, the scattering calculated by Lanczos resampling with closely resembles the sub-box calculation. The Lanczos approach is, however, faster than the method of sub-boxes by approximately three orders of magnitude. This large improvement arises first from the use of the FFT, and second because the scattering is only calculated at the supercell Bragg positions before being resampled to the pixel grid. As the value of is increased, the cutoff function increasingly resembles the rectangular function; hence, by tuning the value of , it is possible to trade off the level of noise that can be tolerated against the amount of blur that is introduced. In particular, choosing ensures that the cutoff function is essentially flat for small , so that correlations at near-neighbour distances can be accurately reflected in the scattering pattern.

I mention three possible modifications to the Lanczos resampling approach. First, if is calculated on a grid with axes parallel to the reciprocal-lattice vectors, equations (45) and (46) can be calculated along each axis in sequence, which reduces the computational expense associated with resampling Getreuer_2011 (). Second, whereas the equations above assume that the idealized cutoff function decays to zero at , it may be preferable to choose a cutoff function that decays to zero at if further noise suppression is required. This corresponds to choosing in the method of sub-boxes. In the Lanczos approach, an idealized cutoff function with that decays to zero at can be obtained for integer by replacing with in equation (46). Finally, it may sometimes be preferable to use a spherically-symmetric cutoff function. This can be achieved by replacing equation (13) with

 (47)

where , is the radius of the largest sphere that can be inscribed in the supercell, and is the supercell volume Marschner_1994 ().

## Vi Conclusions

The main result of this article is that diffuse-scattering calculations from atomistic models can be accelerated by at least two orders of magnitude by applying two fundamental properties of scattering patterns. The first is that the average lattice of a crystal has the periodicity of the crystallographic unit cell, even in the presence of disorder; this property allows the fast Fourier transform (FFT) algorithm to be used to calculate diffuse scattering, which has not previously been considered feasible Butler_1992 (). The second is that the scattering intensity of a periodic supercell is determined at all wavevectors by its values at the Bragg positions of the supercell; this observation allows the scattering intensity to be determined on an arbitrarily fine wavevector grid without redundant calculations. These observations are by no means new: indeed, they date to the earliest work on structural disorder Borie_1971 () and sampling theory in crystallography Sayre_1952 (). Nevertheless, their importance for diffuse-scattering calculations has not previously been recognized.

The acceleration of diffuse-scattering calculations enabled by these results and the associated computer program Scatty (available from the author) has important practical implications. Comparisons of model calculations with the large 3D datasets obtained by modern X-ray and neutron instruments would previously have required many hours, but can now be performed in a few minutes on a laptop. Moreover, it is now practical to fit models of the interactions between disorder variables (local displacements, occupancies, or magnetic moments) directly to 3D experimental datasets in a matter of hours; such fits would previously have required many weeks or months, which has proved prohibitive Welberry_1998 (); Welberry_2008 (). It is also possible to parallelize the intensity calculations for many supercells, potentially leading to further reductions in computer time Gutmann_2010 (). I therefore anticipate that these developments will make quantitative analysis of large diffuse-scattering datasets more practical, more quantitative, and more easily automated. I hope that, similar to the advent of 3D data analysis in conventional crystallography approximately 60 years ago, the advent of 3D diffuse-scattering analysis may allow this approach to be applied to a much wider range of disordered crystalline materials than has yet been possible Welberry_2008 ().

###### Acknowledgements.
I would like to thank Andrew Goodwin and Ross Stewart for valuable discussions related to this work.

### References

1. B. D. Butler, T. R. Welberry, J. Appl. Cryst. 25, 391 (1992).
2. D. A. Keen, J. Phys.: Condens. Matter 14, R819 (2002).
3. T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl, D. Prabhakaran, A. T. Boothroyd, R. J. Aldus, D. F. McMorrow, S. T. Bramwell, Science 326, 415 (2009).
4. J. A. M. Paddison, J. R. Stewart, P. Manuel, P. Courtois, G. J. McIntyre, B. D. Rainford, A. L. Goodwin, Phys. Rev. Lett. 110, 267207 (2013).
5. G. Xu, Z. Zhong, H. Hiraka, G. Shirane, Phys. Rev. B 70, 174109 (2004).
6. M. Paściak, T. R. Welberry, J. Kulda, M. Kempa, J. Hlinka, Phys. Rev. B 85, 224109 (2012).
7. S. Shimomura, N. Wakabayashi, H. Kuwahara, Y. Tokura, Phys. Rev. Lett. 83, 4389 (1999).
8. S. J. L. Billinge, I. Levin, Science 316, 561 (2007).
9. D. A. Keen, A. L. Goodwin, Nature 521, 303 EP (2015).
10. A. R. Overy, A. B. Cairns, M. J. Cliffe, A. Simonov, M. G. Tucker, A. L. Goodwin, Nat. Commun. 7, 10445 EP (2016).
11. T. Weber, A. Simonov, Z. Kristallogr. 227, 238 (2012).
12. D. A. Keen, M. J. Gutmann, C. C. Wilson, J. Appl. Cryst. 39, 714 (2006).
13. F. Ye, Y. Liu, R. Whitfield, R. Osborn, S. Rosenkranz, J. Appl. Cryst. 51, 315 (2018).
14. P. Kraft, A. Bergamaschi, C. Broennimann, R. Dinapoli, E. F. Eikenberry, B. Henrich, I. Johnson, A. Mozzanica, C. M. Schlepütz, P. R. Willmott, B. Schmitt, J. Synchrotron Rad. 16, 368 (2009).
15. M. J. Gutmann, K. Refson, M. v Zimmermann, I. P. Swainson, A. Dabkowski, H. Dabkowska, J. Phys.: Condens. Matter 25, 315402 (2013).
16. T. R. Welberry, T. Proffen, M. Bown, Acta Cryst. A 54, 661 (1998).
17. T. R. Welberry, D. J. Goossens, Acta Cryst. A 64, 23 (2008).
18. T. Weber, H.-B. Bürgi, Acta Cryst. A 58, 526 (2002).
19. R. L. McGreevy, L. Pusztai, Mol. Simul. 1, 359 (1988).
20. M. G. Tucker, D. A. Keen, M. T. Dove, A. L. Goodwin, Q. Hui, J. Phys.: Condens. Matter 19, 335218 (2007).
21. J. A. M. Paddison, A. L. Goodwin, Phys. Rev. Lett. 108, 017204 (2012).
22. J. W. Cooley, J. W. Tukey, Math. Comput. 19, 297 (1965).
23. R. B. Neder, U. Wildgruber, Z. Kristallogr. 186 (1989).
24. S. H. Rahman, Z. Kristallogr. 186, 116 (1989).
25. T. Proffen, R. B. Neder, J. Appl. Cryst. 30, 171 (1997).
26. M. J. Gutmann, J. Appl. Cryst. 43, 250 (2010).
27. F. Frey, Z. Kristallogr. 212, 257 (1997).
28. E. Whittaker, P. Roy. Soc. Edinb. A 35, 181 (1915).
29. C. E. Shannon, Proc. IRE 37, 10 (1949).
30. D. Sayre, Acta Cryst. 5, 843 (1952).
31. R. Singleton, IEEE Trans. Audio and Electroacoustics 17, 93 (1969).
32. D. J. Goossens, A. P. Heerdegen, E. J. Chan, T. R. Welberry, Metal. Mater. Trans. A 42, 23 (2011).
33. T. M. Michels-Clark, Methods for quantitative local structure analysis of crystalline materials employing high performance computing, Ph.D. thesis, University of Tennessee–Knoxville (2014).
34. R. Withers, IUCrJ 2, 74 (2015).
35. B. Borie, C. J. Sparks, Jnr, Acta Cryst. A 27, 198 (1971).
36. M. Hayakawa, J. B. Cohen, Acta Cryst. A 31, 635 (1975).
37. G. L. Squires, Introduction to the Theory of Thermal Neutron Scattering (Cambridge University Press, Cambridge, 1978).
38. S. W. Lovesey, Theory of Neutron Scattering from Condensed Matter: Polarization Effects and Magnetic Scattering, vol. 2 (Oxford University Press, Oxford, 1987).
39. P. H. Conlon, Aspects of frustrated magnetism, Ph.D. thesis, University of Oxford (2010).
40. S. T. Bramwell, M. J. Gingras, Science 294, 1495 (2001).
41. L. Pauling, J. Am. Chem. Soc. 57, 2680 (1935).
42. W. F. Giauque, J. W. Stout, J. Am. Chem. Soc. 58, 1144 (1936).
43. A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, B. S. Shastry, Nature 399, 333 (1999).
44. J. S. Gardner, B. D. Gaulin, S.-H. Lee, C. Broholm, N. P. Raju, J. E. Greedan, Phys. Rev. Lett. 83, 211 (1999).
45. P. M. M. Thygesen, J. A. M. Paddison, R. Zhang, K. A. Beyer, K. W. Chapman, H. Y. Playford, M. G. Tucker, D. A. Keen, M. A. Hayward, A. L. Goodwin, Phys. Rev. Lett. 118, 067201 (2017).
46. B. Wehinger, D. Chernyshov, M. Krisch, S. Bulat, V. Ezhov, A. Bosak, J. Phys.: Condens. Matter 26, 265401 (2014).
47. B. D. Butler, T. R. Welberry, Acta Cryst. A 49, 736 (1993).
48. J. A. M. Paddison, J. R. Stewart, A. L. Goodwin, J. Phys.: Condens. Matter 25, 454220 (2013).
49. J. W. Gibbs, Nature 59, 606 (1899).
50. C. Lanczos, Applied Analysis (Prentice-Hall, Englewood Cliffs, New Jersey, 1956).
51. C. E. Duchon, J. Appl. Meteorol. 18, 1016 (1979).
52. P. Getreuer, IPOL Image Proc. On Line 1, 238 (2011).
53. S. R. Marschner, R. J. Lobb, Proceedings of the Conference on Visualization ’94, VIS ’94 (IEEE Computer Society Press, Los Alamitos, CA, USA, 1994), pp. 100–107.
54. J. A. M. Paddison, M. J. Gutmann, J. R. Stewart, M. G. Tucker, M. T. Dove, D. A. Keen, A. L. Goodwin, Phys. Rev. B 97, 014429 (2018).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters