Integration, scaling, space-group assignment and post refinement

Kabsch, W.

doi:10.1107/97809553602060000676

International
Tables for
Crystallography
Volume F
Crystallography of biological macromolecules
Edited by M. G. Rossmann and E. Arnold

pdf | chapter contents | chapter index | related articles

International Tables for Crystallography (2006). Vol. F. ch. 11.3, pp. 218-225 | 1 | 2 |
https://doi.org/10.1107/97809553602060000676

Chapter 11.3. Integration, scaling, space-group assignment and post refinement

W. Kabsch^a ^*

^a Max-Planck-Institut für medizinische Forschung, Abteilung Biophysik, Jahnstrasse 29, 69120 Heidelberg, Germany
Correspondence e-mail: kabsch@mpimf-heidelberg.mpg.de

The key steps in the processing of diffraction data from single crystals are described. The topics covered include: the modelling of the positions of all the reflections recorded in the images; the integration of diffraction intensities; data correction, scaling and post refinement; and space-group assignment. The principles of the methods are described as they are employed by the program XDS (Section 25.2.9 ).

Keywords: Bravais lattice; background; basis vectors; data integration; data processing; diffraction-pattern prediction; indexing; integration of diffraction data; partiality; post refinement; profile fitting; reflection centroid; scaling; space-group assignment; standard profiles.

11.3.1. Introduction

| top | pdf |

Key steps in the processing of diffraction data from single crystals involve: (a) accurate modelling of the positions of all the reflections recorded in the images; (b) integration of diffraction intensities; (c) data correction, scaling and post refinement; and (d) space-group assignment. Much of the theory and many of the methods for carrying out these steps were developed about two decades ago for processing rotation data recorded on film and were later extended to exploit fully the capabilities of a variety of electronic area detectors; some CCD (charge-coupled device) and multiwire detectors allow the recording of finely sliced rotation data because of their fast data read-out. In this chapter, the principles of the methods are described as they are employed by the program XDS (Section 25.2.9 ). These apply equally well to rotation images covering small or large oscillation ranges. A large number of other systems have been developed which differ in the details of the implementations. Some of these packages are described in Chapter 25.2 . The theory and practice of processing fine-sliced data have recently been discussed by Pflugrath (1997).

11.3.2. Modelling rotation images

| top | pdf |

The observed diffraction pattern, i.e., the positions of the reflections recorded in the rotation-data images, is controlled by a small set of parameters which must be accurately determined before integration can start. Approximate values for some of these parameters are given by the experimental setup, whereas others may be completely unknown and must be obtained from the rotation images. This is achieved by automatic location of strong diffraction spots, extraction of a primitive lattice basis that yields integer indices for the observed reflections, and subsequent refinement of all parameters to minimize the discrepancies between observed and calculated spot positions in the data images.

11.3.2.1. Coordinate systems and parameters

| top | pdf |

In the rotation method, the incident beam wave vector $[{\bf S}_{0}]$ of length $[1/\lambda]$ (λ is the wavelength) is fixed while the crystal is rotated around a fixed axis described by a unit vector $[{\bf m}_{2}]$ . $[{\bf S}_{0}]$ points from the X-ray source towards the crystal. It is assumed that the incident beam and the rotation axis intersect at one point at which the crystal must be located. This point is defined as the origin of a right-handed orthonormal laboratory coordinate system $[\{{\bf l}_{1},{\bf l}_{2},{\bf l}_{3}\}]$ . This fixed but otherwise arbitrary system is used as a reference frame to specify the setup of the diffraction experiment.

Diffraction data are assumed to be recorded on a fixed planar detector. A right-handed orthonormal detector coordinate system $[\{{\bf d}_{1},{\bf d}_{2},{\bf d}_{3}\}]$ is defined such that a point with coordinates X, Y in the detector plane is represented by the vector $[(X - X_{0}){\bf d}_{1} + {(Y - Y_{0}){\bf d}_{2} + F{\bf d}_{3}}]$ with respect to the laboratory coordinate system. The origin $[X_{0}, Y_{0}]$ of the detector plane is found at a distance [|F|] from the crystal position. It is assumed that the diffraction data are recorded on adjacent non-overlapping rotation images, each covering a constant oscillation range $[\Delta_{\varphi}]$ with image No. 1 starting at spindle angle $[\varphi_{0}]$ .

Diffraction geometry is conveniently expressed with respect to a right-handed orthonormal goniostat system $[\{{\bf m}_{1}, {\bf m}_{2}, {\bf m}_{3}\}]$ . It is constructed from the rotation axis and the incident beam direction such that $[{\bf m}_{1} = ({\bf m}_{2}\times {\bf S}_{0})/|{\bf m}_{2} \times {\bf S}_{0}|]$ and $[{\bf m}_{3} = {\bf m}_{1} \times {\bf m}_{2}]$ . The origin of the goniostat system is defined to coincide with the origin of the laboratory system.

Finally, a right-handed crystal coordinate system $[\{{\bf b}_{1}, {\bf b}_{2}, {\bf b}_{3}\}]$ and its reciprocal basis $[\{{\bf b}_{1}^{*}, {\bf b}_{2}^{*}, {\bf b}_{3}^{*}\}]$ are defined to represent the unrotated crystal, i.e., at rotation angle $[\varphi = 0^{\circ}]$ , such that any reciprocal-lattice vector can be expressed as $[{\bf p}_{0}^{*} = h{\bf b}_{1}^{*} + k{\bf b}_{2}^{*} + l{\bf b}_{3}^{*}]$ where [h, k, l] are integers.

Using a Gaussian model, the shape of the diffraction spots is specified by two parameters: the standard deviations of the reflecting range $[\sigma_{M}]$ and the beam divergence $[\sigma_{D}]$ (see Section 11.3.2.3 ). This leads to an integration region around the spot defined by the parameters $[\delta_{M}]$ and $[\delta_{D}]$ , which are typically chosen to be 6–10 times larger than $[\sigma_{M}]$ and $[\sigma_{D}]$ , respectively.

Knowledge of the parameters $[{\bf S}_{0}]$ , $[{\bf m}_{2}]$ , $[{\bf b}_{1}]$ , $[{\bf b}_{2}]$ , $[{\bf b}_{3}]$ , $[X_{0}]$ , $[Y_{0}, F]$ , $[{\bf d}_{1}]$ , $[{\bf d}_{2}]$ , $[{\bf d}_{3}]$ , $[\varphi_{0}]$ and $[\Delta_{\varphi}]$ is sufficient to compute the location of all diffraction peaks recorded in the data images. Determination and refinement of these parameters are described in the following sections.

11.3.2.2. Spot prediction

| top | pdf |

It is assumed here that accurate values of all parameters describing the diffraction experiment are available, permitting prediction of the positions of all diffraction peaks recorded in the data images. Let $[{\bf p}^{*}_{0}]$ denote any arbitrary reciprocal-lattice vector if the crystal has not been rotated, i.e., at rotation angle $[\varphi = 0^{\circ}]$ . $[{\bf p}^{*}_{0}]$ can be expressed by its components with respect to the orthonormal goniostat system as $[{\bf p}^{*}_{0} = {\bf m}_{1} ({\bf m}_{1} \cdot {\bf p}^{*}_{0}) + {\bf m}_{2}({\bf m}_{2} \cdot {\bf p}^{*}_{0}) + {\bf m}_{3}({\bf m}_{3} \cdot {\bf p}^{*}_{0}).]$ Depending on the diffraction geometry, $[{\bf p}^{*}_{0}]$ may be rotated into a position fulfilling the reflecting condition. The required rotation angle φ and the coordinates X, Y of the diffracted beam at its intersection with the detector plane can be found from $[{\bf p}^{*}_{0}]$ as follows.

Rotation by φ around axis $[{\bf m}_{2}]$ changes $[{\bf p}^{*}_{0}]$ into $[{\bf p}^{*}]$ . $[\eqalign{{\bf p}^{*} &=D({\bf m}_{2}, \varphi){\bf p}^{*}_{0} = {\bf m}_{2}({\bf m}_{2} \cdot {\bf p}^{*}_{0}) + [{\bf p}^{*}_{0} - {\bf m}_{2}({\bf m}_{2} \cdot {\bf p}^{*}_{0})]\cos \varphi\cr&\quad+ \;{\bf m}_{2} \times {\bf p}^{*}_{0} \sin \varphi\cr &= {\bf m}_{1}({\bf m}_{1} \cdot {\bf p}^{*}_{0} \cos \varphi + {\bf m}_{3} \cdot {\bf p}^{*}_{0} \sin \varphi) + {\bf m}_{2}{\bf m}_{2} \cdot {\bf p}^{*}_{0}\cr &\quad+\; {\bf m}_{3}({\bf m}_{3} \cdot {\bf p}^{*}_{0}\cos \varphi - {\bf m}_{1} \cdot {\bf p}^{*}_{0}\sin \varphi)\cr &= {\bf m}_{1}({\bf m}_{1} \cdot {\bf p}^{*}) + {\bf m}_{2}({\bf m}_{2} \cdot {\bf p}^{*}) + {\bf m}_{3}({\bf m}_{3} \cdot {\bf p}^{*}).}]$ The incident and diffracted beam wave vectors, $[{\bf S}_{0}]$ and S, have their termini on the Ewald sphere and satisfy the Laue equations $[{\bf S} = {\bf S}_{0} + {\bf p}^{*}, \quad {\bf S}^{2} = {\bf S}_{0}^{2}\;\Longrightarrow\; {\bf p}^{*2} = -2{\bf S}_{0} \cdot {\bf p}^{*} = {\bf p}^{*2}_{0}.]$ If $[\rho = [{\bf p}^{*2}_{0} - ({\bf p}^{*}_{0} \cdot {\bf m}_{2})^{2}]^{1/2}]$ denotes the distance of $[{\bf p}^{*}_{0}]$ from the rotation axis, solutions for $[{\bf p}^{*}]$ and φ can be obtained in terms of $[{\bf p}^{*}_{0}]$ as $[\eqalign{{\bf p}^{*} \cdot {\bf m}_{3} &= [-{\bf p}^{*2}_{0}/2 - ({\bf p}^{*}_{0} \cdot {\bf m}_{2})({\bf S}_{0} \cdot {\bf m}_{2})]/{\bf S}_{0}\cdot {\bf m}_{3}\cr {\bf p}^{*} \cdot {\bf m}_{2} &= {\bf p}^{*}_{0} \cdot {\bf m}_{2}\cr {\bf p}^{*} \cdot {\bf m}_{1} &= \pm [\rho^{2} - ({\bf p}^{*} \cdot {\bf m}_{3})^{2}]^{1/2}\cr \cos \varphi &= [({\bf p}^{*} \cdot {\bf m}_{1})({\bf p}^{*}_{0} \cdot {\bf m}_{1}) + ({\bf p}^{*} \cdot {\bf m}_{3})({\bf p}^{*}_{0} \cdot {\bf m}_{3})]/\rho^{2}\cr \sin \varphi &= [({\bf p}^{*} \cdot {\bf m}_{1})({\bf p}^{*}_{0} \cdot {\bf m}_{3}) - ({\bf p}^{*} \cdot {\bf m}_{3})({\bf p}^{*}_{0} \cdot {\bf m}_{1})]/\rho^{2}.}]$ In general, there are two solutions according to the sign of $[{\bf p}^{*} \cdot {\bf m}_{1}]$ . If $[\rho^{2} \lt ({\bf p}^{*} \cdot {\bf m}_{3})^{2}]$ or $[{\bf p}^{*2}_{0}\gt 4{\bf S}_{0}^{2}]$ , the Laue equations have no solution and the reciprocal-lattice point $[{\bf p}^{*}_{0}]$ is in the `blind' region.

If $[F{\bf S} \cdot {\bf d}_{3} \gt 0]$ , the diffracted beam intersects the detector plane at the point $[\eqalign{F{\bf S}/{\bf S} \cdot {\bf d}_{3} &= (F{\bf S} \cdot {\bf d}_{1}/{\bf S} \cdot {\bf d}_{3}){\bf d}_{1} + (F{\bf S} \cdot {\bf d}_{2}/{\bf S} \cdot {\bf d}_{3}){\bf d}_{2} + F{\bf d}_{3}\cr &= (X - X_{0}){\bf d}_{1} + (Y - Y_{0}){\bf d}_{2} + F{\bf d}_{3},}]$ which leads to a diffraction spot recorded at detector coordinates $[\eqalign{X &= X_{0} + F{\bf S} \cdot {\bf d}_{1}/{\bf S} \cdot {\bf d}_{3},\cr Y &= Y_{0} + F{\bf S} \cdot {\bf d}_{2}/{\bf S} \cdot {\bf d}_{3}.}]$

11.3.2.3. Standard spot shape

| top | pdf |

A reciprocal-lattice point crosses the Ewald sphere by the shortest route only if the crystal happens to be rotated about an axis perpendicular to both the diffracted and incident beam wave vectors, the `β-axis' $[{\bf e}_{1} = {\bf S} \times {\bf S}_{0}/|{\bf S} \times {\bf S}_{0}|]$ , as introduced by Schutt & Winkler (1977 ). Rotation around the fixed axis $[{\bf m}_{2}]$ , as enforced by the rotation camera, thus leads to an increase in the length of the shortest path by the factor $[1/|{\bf e}_{1} \cdot {\bf m}_{2}|]$ . This has motivated the introduction of a coordinate system $[\{{\bf e}_{1}, {\bf e}_{2}, {\bf e}_{3}\}]$ , specific for each reflection, which has its origin on the surface of the Ewald sphere at the terminus of the diffracted beam wave vector S, $[\displaylines{{\bf e}_{1} = {\bf S} \times {\bf S}_{0}/|{\bf S} \times {\bf S}_{0}|, \quad {\bf e}_{2} = {\bf S} \times {\bf e}_{1}/|{\bf S} \times {\bf e}_{1}|,\cr \quad {\bf e}_{3} = ({\bf S} + {\bf S}_{0})/|{\bf S} + {\bf S}_{0}|.}]$ The unit vectors $[{\bf e}_{1}]$ and $[{\bf e}_{2}]$ are tangential to the Ewald sphere, while $[{\bf e}_{3}]$ is perpendicular to $[{\bf e}_{1}]$ and $[{\bf p}^{*} = {\bf S} - {\bf S}_{0}]$ . The shape of a reflection, as represented with respect to $[\{{\bf e}_{1}, {\bf e}_{2}, {\bf e}_{3}\}]$ , then no longer contains geometrical distortions resulting from the fixed rotation axis of the camera and the oblique incidence of the diffracted beam on a flat detector. Instead, all reflections appear as if they had followed the shortest path through the Ewald sphere and had been recorded on the surface of the sphere.

A detector pixel at X′, Y′ in the neighbourhood of the reflection centre X, Y, when the crystal is rotated by φ′ instead of φ, is mapped to the profile coordinates $[\varepsilon_{1}, \varepsilon_{2}, \varepsilon_{3}]$ by the following procedure: $[\eqalign{{\bf S}' &= [(X' - X_{0}){\bf d}_{1} + (Y' - Y_{0}){\bf d}_{2} + F{\bf d}_{3}]\cr &\quad\times \{\lambda \cdot [(X' - X_{0})^{2} + (Y' - Y_{0})^{2} + F^{2}]^{1/2}\}^{-1}\cr \varepsilon_{1} &= {\bf e}_{1} \cdot ({\bf S}' - {\bf S})180/(|{\bf S}|\pi ),\cr \varepsilon_{2} &= {\bf e}_{2} \cdot ({\bf S}' - {\bf S})180/(|{\bf S}|\pi )\cr \varepsilon_{3} &= {\bf e}_{3} \cdot [D({\bf m}_{2}, \varphi' - \varphi){\bf p}^{*} - {\bf p}^{*}]180/(|{\bf p}^{*}|\pi ) \simeq \zeta \cdot (\varphi' - \varphi)\cr \zeta &= {\bf m}_{2} \cdot {\bf e}_{1}.}]$ ζ corrects for the increased path length of the reflection through the Ewald sphere and is closely related to the reciprocal Lorentz correction factor $[L^{-1} = |{\bf m}_{2} \cdot ({\bf S} \times {\bf S}_{0})|/(|{\bf S}| \cdot |{\bf S}_{0}|) = |\zeta \cdot \sin \angle ({\bf S},{\bf S}_{0})|.]$

Because of crystal mosaicity and beam divergence, the intensity of a reflection is smeared around the diffraction maximum. The fraction of total reflection intensity found in the volume element $[\hbox{d}\varepsilon_{1} \hbox{d}\varepsilon_{2} \hbox{d}\varepsilon_{3}]$ at $[\varepsilon_{1}, \varepsilon_{2}, \varepsilon_{3}]$ can be approximated by Gaussian functions: $[\displaylines{\omega(\varepsilon_{1}, \varepsilon_{2}, \varepsilon_{3})\hbox{d}\varepsilon_{1} \hbox{d}\varepsilon_{2} \hbox{d}\varepsilon_{3}\hfill\cr \quad= {\exp(- \varepsilon_{1}^{2}/2\sigma_{D}^{2}) \over (2\pi )^{1/2} \sigma_{D}} \hbox{d}\varepsilon_{1} \cdot {\exp(- \varepsilon_{2}^{2}/2\sigma_{D}^{2}) \over (2\pi )^{1/2} \sigma_{D}} \hbox{d}\varepsilon_{2} \cdot {\exp(- \varepsilon_{3}^{2}/2\sigma_{M}^{2}) \over (2\pi )^{1/2} \sigma_{M}} \hbox{d}\varepsilon_{3}.}]$

11.3.2.4. Spot centroids and partiality

| top | pdf |

The intensity of a reflection can be completely recorded on one image, or distributed among several adjacent images. The fraction $[R_{j}]$ of total intensity recorded on image j, the `partiality' of the reflection, can be derived from the distribution function $[\omega (\varepsilon_{1}, \varepsilon_{2}, \varepsilon_{3})]$ as $[\eqalign{R_{j} &= {\textstyle\int\limits_{- \infty}^{\infty}} \hbox{d}\varepsilon_{1} {\textstyle\int\limits_{- \infty}^{\infty}} \hbox{d}\varepsilon_{2} {\textstyle\int\limits_{\zeta [\varphi_{0} + (j - 1)\Delta_{\varphi} - \varphi]}^{\zeta (\varphi_{0} + j\Delta_{\varphi} - \varphi)}} \hbox{d}\varepsilon_{3} \;\omega(\varepsilon_{1}, \varepsilon_{2}, \varepsilon_{3})\cr &= \{[1/(2\pi )^{1/2} \sigma_{M}]/|\zeta|\} \cr&\quad\times{\textstyle\int\limits_{\varphi_{0} + (j - 1)\Delta_{\varphi}}^{\varphi_{0} + j\Delta_{\varphi}}} \exp[- (\varphi' - \varphi)^{2}/2(\sigma_{M}/|\zeta|)^{2}] \;\hbox{d}\varphi'\cr &= \left(\hbox{erf}[|\zeta|(\varphi_{0} + j\Delta_{\varphi} - \varphi)/(2)^{1/2} \sigma_{M}]\right. \cr&\left.\quad-\;\hbox{erf}\{|\zeta|[\varphi_{0} + (j - 1)\Delta_{\varphi} - \varphi]/(2)^{1/2} \sigma_{M}\}\right)\big/2.}]$ The integral is evaluated by using a numerical approximation of the error function, erf (Abramowitz & Stegun, 1972 ).

While the spot centroids in the detector plane are usually good estimates for the detector position of the diffraction maximum, the angular centroid about the rotation axis, $[Z = \varphi_{0} + \Delta_{\varphi} \cdot {\textstyle\sum\limits_{j = -\infty}^{\infty}} (j - 1/2) R_{j} \approx \varphi,]$ can be a rather poor guess for the true φ angle of the maximum. Its accuracy depends strongly on the value of φ and the size of the oscillation range $[\Delta_{\varphi}]$ relative to the mosaicity $[\sigma_{M}]$ of the crystal. For a reflection fully recorded on image j, the value $[Z = \varphi_{0} + {(\;j - 1/2) \cdot \Delta_{\varphi}}]$ will always be obtained, which is correct only if φ accidentally happens to be close to the centre of the rotation range of the image. In contrast, the φ angle of a partial reflection recorded on images j and [j + 1] is closely approximated by $[Z = \varphi_{0} + [\;j + (R_{j + 1} - R_{j})/2] \cdot \Delta_{\varphi}]$ . If many images contribute to the spot intensity, $[Z(\varphi)]$ is always an excellent approximation to the ideal angular position φ when the Laue equations are satisfied; in fact, in the limiting case of infinitely fine-sliced data, it can be shown that $[\lim_{\Delta_{\varphi} \rightarrow 0}Z(\varphi) = \varphi]$ .

Most refinement routines minimize the discrepancies between the predicted φ angles and their approximations obtained from the observed Z centroids, and must therefore carefully distinguish between fully and partially recorded reflections. This distinction is unnecessary, however, if observed Z centroids are compared with their analytic forms instead, because the sensitivity of the centroid positions to the diffraction parameters is correctly weighted in either case (see Section 11.3.2.8 ).

11.3.2.5. Localizing diffraction spots

| top | pdf |

Recognition and refinement of the parameter values controlling the observed diffraction pattern begins with the extraction of a list of coordinates of strong spots occurring in the images. As implemented in XDS, this list is obtained by the following procedure. First, each pixel value is compared with the mean value and standard deviation of surrounding pixels in the same image and classified as a strong pixel if its value exceeds the mean by a given multiple (typically 3 to 5) of the standard deviation. Values of the strong pixels and their location addresses and image running numbers are stored in a hash table during spot search [for a discussion of the hash technique, see Wirth (1976 )]. After processing a fixed number of images, or when the table is full, all strong pixels are labelled by a unique number identifying the spot to which they belong. By definition, any two such pixels which can be connected by direct strong neighbours in two or three dimensions (if there are adjacent images) belong to the same spot (equivalence class). The labelling is achieved by the highly efficient algorithm for the recording of equivalence classes developed by Rem (see Dijkstra, 1976 ). At the end of this procedure, the table is searched for spots that have no contributing strong pixel on the current or the previous image. These spots are complete and their centroids are evaluated and saved in a file. To make room for new strong pixels as the spot search proceeds, all entries of strong pixels that are no longer needed are removed from the hash table and the remaining ones are rehashed. On termination, a list $[X'_{i}, Y'_{i}, Z'_{i}\ (i = 1, \ldots, n)]$ of the centroids of strong spots is available.

11.3.2.6. Basis extraction

| top | pdf |

Any reciprocal-lattice vector can be written in the form $[{\bf p}_{0}^{*} = h{\bf b}_{1}^{*} + k{\bf b}_{2}^{*} + l{\bf b}_{3}^{*}]$ where h, k, l are integer numbers and $[{\bf b}_{1}^{*}, {\bf b}_{2}^{*}, {\bf b}_{3}^{*}]$ are basis vectors of the lattice. The basis vectors which describe the orientation, metric and symmetry of the crystal, as well as the reflection indices h, k, l, have to be determined from the list of strong diffraction spots $[X'_{i}, Y'_{i}, Z'_{i}\ (i = 1, \ldots, n)]$ . Ideally, each spot corresponds to a reciprocal-lattice vector $[{\bf p}_{0}^{*}]$ which satisfies the Laue equations after a crystal rotation by φ. Substituting the observed value Z′ for the unknown φ angle (see Section 11.3.2.4 ), $[{\bf p}_{0}^{*}]$ is found from the observed spot coordinates as $[\eqalign{{\bf p}_{0}^{*} &=D({\bf m}_{2}, -Z') ({\bf S'} - {\bf S}_{0})\cr {\bf S'} &= \left[(X'- X_{0}){\bf d}_{1} + (Y'- Y_{0}){\bf d}_{2} + F{\bf d}_{3}\right]\cr &\quad\times\left\{\lambda \cdot \left[(X'- X_{0})^{2} + (Y' - Y_{0})^{2} + F^{2}\right]^{1/2}\right\}^{-1}.}]$ Unfortunately, the reciprocal-lattice vectors $[{\bf p}^{*}_{0i}\ (i = 1, \ldots, n)]$ derived from the above list of strong diffraction spots often contain a number of `aliens' (spots arising from fluctuations of the background, from ice, or from satellite crystals) and a robust method has to be used which is still capable of recognizing the dominant lattice. One approach, suggested by Bricogne (1986 ) and implemented in a number of variants (Otwinowski & Minor, 1997 ; Steller et al., 1997 ), is to identify a lattice basis as the three shortest linear independent vectors $[{\bf b}_{1}, {\bf b}_{2}, {\bf b}_{3}]$ , each at a maximum of the Fourier transform $[\sum\nolimits_{i=1}^{n} \cos (2\pi {\bf b} \cdot {\bf p}^{*}_{0i})]$ . Alternatively, a reciprocal basis for the dominant lattice can be determined from short differences between the reciprocal-lattice vectors (Howard, 1986 ; Kabsch, 1988a ). As implemented in XDS, a lattice basis is found by the following procedure.

The list of given reciprocal-lattice points $[{\bf p}^{*}_{0i}\ (i = 1, \ldots,n)]$ is first reduced to a small number m of low-resolution difference-vector clusters $[{\bf v}^{*}_{\mu}\ (\mu = 1, \ldots, m)]$ . $[f_{\mu}]$ is the population of a difference-vector cluster $[{\bf v}^{*}_{\mu}]$ , that is the number of times the difference between any two reciprocal-lattice vectors $[{\bf p}^{*}_{0i} - {\bf p}^{*}_{0j}]$ is approximately equal to $[{\bf v}^{*}_{\mu}]$ . In a second step, three linear independent vectors $[{\bf b}_{1}^{*}, {\bf b}_{2}^{*}, {\bf b}_{3}^{*}]$ are selected among all possible triplets of difference-vector clusters that maximize the function Q: $[\eqalign{Q({\bf b}_{1}^{*}, {\bf b}_{2}^{*}, {\bf b}_{3}^{*}) &={\textstyle\sum\limits_{\mu = 1}^{m}} f_{\mu} q(\xi_{1}^{\mu}, \xi_{2}^{\mu}, \xi_{3}^{\mu})\cr q(\xi_{1}^{\mu}, \xi_{2}^{\mu}, \xi_{3}^{\mu}) &= \exp\left(-2 {\textstyle\sum\limits_{k = 1}^{3}}\left\{[\max(|\xi_{k}^{\mu} - h_{k}^{\mu}| - \varepsilon,0)/\varepsilon]^{2}\right.\right.\cr&\quad\left. +\; [\max(|h_{k}^{\mu}| - \delta,0)]^{2}\right\}\bigg)\cr \xi_{k}^{\mu} &={\bf v}_{\mu}^{*} \cdot {\bf b}_{k}, \quad {\bf v}_{\mu}^{*} = {\textstyle\sum\limits_{k = 1}^{3}} \xi_{k}^{\mu}{\bf b}_{k}^{*}, \quad {\bf b}_{k} \cdot {\bf b}_{l}^{*} = {\left\{\matrix{1 \hbox{ if } k = l \hbox{;}\hfill\cr 0 \hbox{ otherwise}\cr}\right.}\cr h_{k}^{\mu} &= \hbox{ nearest integer to}\ \xi_{k}^{\mu}.}\hfill]$ The absolute maximum of Q is assumed if all difference vectors can be expressed as small integral multiples of the best triplet. Deviations from this ideal situation are quantified by the quality measure q. The value of q declines sharply if the expansion coefficients $[\xi_{k}^{\mu}]$ deviate by more than ɛ from their nearest integers $[h_{k}^{\mu}]$ or if the indices are absolutely larger than δ. The constraint on the allowed range of indices prevents the selection of a spurious triplet of very short difference vector clusters which might be present in the set. Excellent results have been obtained using $[\varepsilon = 0.05]$ and $[\delta = 5]$ . The best vector triplet thus found is refined against the observed difference-vector clusters. Finally, a reduced cell is derived from the refined reciprocal-base vector triplet as defined in IT A (2005 ), Chapter 9.2 .

11.3.2.7. Indexing

| top | pdf |

Once a basis $[{\bf b}_{1}, {\bf b}_{2}, {\bf b}_{3}]$ of the lattice is available, integral indices $[h_{i}, k_{i}, l_{i}]$ must be assigned to each reciprocal-lattice vector $[{\bf p}_{0i}^{*}\ (i = 1, \ldots, n)]$ . Using the integers nearest to $[{\bf p}_{0i}^{*} \cdot {\bf b}_{k}\ (k = 1,2,3)]$ as indices of the reciprocal-lattice vectors $[{\bf p}_{0i}^{*}]$ could easily lead to a misindexing of longer vectors because of inaccuracies in the basis vectors $[{\bf b}_{k}]$ and the initial values of the parameters describing the instrumental setup. A more robust solution of the indexing problem is provided by the local indexing method which assigns only small index differences $[h_{i} - h_{j}, k_{i} - k_{j}, l_{i} - l_{j}]$ between pairs of neighbouring reciprocal-lattice vectors (Kabsch, 1993 ).

The reciprocal-lattice points can be considered as the nodes of a tree. The tree connects the n points to each other with the connections as its branches. The length $[\ell_{ij}]$ of a possible branch between nodes i and j is defined here as $[\eqalign{\ell_{ij} &= 1 - \exp\left(-2 {\textstyle\sum\limits_{k = 1}^{3}}\left\{[\max(|\xi_{k}^{ij} - h_{k}^{ij}| - \varepsilon,0)/\varepsilon]^{2}\right.\right.\cr&\quad\left. +\; [\max(|h_{k}^{ij}| - \delta,0)]^{2}\right\}\bigg)\cr \xi_{k}^{ij} &= ({\bf p}_{0i}^{*} - {\bf p}_{0j}^{*}) \cdot {\bf b}_{k}, \quad h_{k}^{ij} = \hbox{ nearest integer of }\xi_{k}^{ij}, \quad k = 1,2,3.}]$ Reliable index differences are indicated by short branches; in fact, $[\ell_{ij}]$ is 0 if none of the indices $[h_{k}^{ij}]$ is absolutely larger than δ and the $[\xi_{k}^{ij}]$ are integer values to within ɛ. Typical values of ɛ and δ are $[\varepsilon = 0.05]$ and $[\delta = 5]$ . Defining the length of a tree as the sum of the lengths of its branches, a shortest tree among all $[n^{n-2}]$ possible trees is determined by the elegant algorithm described by Dijkstra (1976 ). Starting with arbitrary indices 0, 0, 0 for the root node, the local indexing method then consists of traversing the shortest tree and thereby assigning each node the indices of its predecessor plus the small index differences between the two nodes.

During traversal of the tree, each node is also given a subtree number. Starting with subtree number 1 for the root node, each successor node is given the same subtree number as its predecessor if the length of the connecting branch is below a minimal length $[\ell_{\min}]$ . Otherwise its subtree number is incremented by 1. Thus all nodes in the same subtree have internally consistent reflection indices. Defining the size of a subtree by the number of its nodes, aliens are usually found in small subtrees. Finally, a constant index offset is determined such that the centroids of the observed reciprocal-lattice points $[{\bf p}_{0i}^{*}]$ belonging to the largest subtree and their corresponding grid vectors $[\sum\nolimits_{k = 1}^{3}h_{k}^{i}{\bf b}_{k}^{*}]$ are as close as possible. This offset is added to the indices of each reciprocal-lattice point.

11.3.2.8. Refinement

| top | pdf |

For a fixed detector, the diffraction pattern depends on the parameters $[{\bf S}_{0}, {\bf m}_{2}, {\bf b}_{1}, {\bf b}_{2}, {\bf b}_{3}, X_{0}, Y_{0}]$ and F. Starting values for the parameters can be obtained by the procedures described above that do not rely on prior knowledge of the crystal orientation, space-group symmetry or unit-cell metric. Better estimates of the parameter values, as required for the subsequent integration step, can be obtained by the method of least squares from the list of n observed indexed reflection centroids $[h_{i}]$ , $[k_{i}]$ , $[l_{i}]$ , $[X_{i}']$ , $[Y_{i}']$ , $[Z_{i}']$ $[(i = 1, \ldots, n)]$ . In this method, the parameters are chosen to minimize a weighted sum of squares of the residuals $[E = w_{X} {\textstyle\sum\limits_{i=1}^{n}} (\Delta_{X}^{i})^{2} + w_{Y} {\textstyle\sum\limits_{i=1}^{n}} (\Delta_{Y}^{i})^{2} + w_{Z} {\textstyle\sum\limits_{i=1}^{n}} (\Delta_{Z}^{i})^{2}.]$ The residuals between the calculated $[(X_{i}, Y_{i}, Z_{i})]$ and observed spot centroids are $[\eqalign{\Delta_{X}^{i} &=X_{i} - X_{i}' =X_{0} + F{\bf S}_{i} \cdot {\bf d}_{1}/{\bf S}_{i} \cdot {\bf d}_{3} - X_{i}'\cr \Delta_{Y}^{i} &=Y_{i} - Y_{i}' =Y_{0} + F{\bf S}_{i} \cdot {\bf d}_{2}/{\bf S}_{i} \cdot {\bf d}_{3} - Y_{i}'\cr \Delta_{Z}^{i} &=Z_{i} - Z_{i}' =\varphi_{0} + \Delta_{\varphi} {\textstyle\sum\limits_{j = -\infty}^{\infty}} (\;j - 1/2) R_{j}^{i} - Z_{i}'.}]$

Let $[s_{\mu}\ (\mu = 1, \ldots, k)]$ denote the k independent parameters for which initial estimates are available. Expanding the residuals to first order in the parameter changes $[\delta s_{\mu}]$ gives $[\Delta (s_{\mu} + \delta s_{\mu}) \approx \Delta (s_{\mu}) + \sum\limits_{\mu = 1}^{k} {\partial\Delta\over \partial s_{\mu}}\delta s_{\mu}.]$ The parameters should be changed in such a way as to minimize $[E(\delta s_{\mu})]$ , which implies $[\partial E/\partial \delta s_{\mu} = 0]$ for $[\mu = 1, \ldots, k]$ . The $[\delta s_{\mu}]$ are found as the solution of the k normal equations $[\eqalign{&\sum\limits_{\mu'= 1}^{k} \left(w_{X} \sum\limits_{i=1}^{n} {\partial\Delta_{X}^{i} \over \partial s_{\mu}} {\partial\Delta_{X}^{i} \over \partial s_{\mu'}} + w_{Y}\sum\limits_{i=1}^{n} {\partial\Delta_{Y}^{i} \over \partial s_{\mu}} {\partial\Delta_{Y}^{i} \over \partial s_{\mu'}} + w_{Z} \sum\limits_{i=1}^{n} {\partial\Delta_{Z}^{i} \over \partial s_{\mu}} {\partial\Delta_{Z}^{i} \over \partial s_{\mu'}}\right) \delta s_{\mu'}\cr &\quad = - \left(w_{X}\sum\limits\limits_{i=1}^{n} \Delta_{X}^{i} {\partial \Delta_{X}^{i} \over \partial s_{\mu}} + w_{Y}\sum\limits\limits_{i=1}^{n}\Delta_{Y}^{i} {\partial \Delta_{Y}^{i} \over \partial s_{\mu}} + w_{Z}\sum\limits\limits_{i=1}^{n} \Delta_{Z}^{i} {\partial \Delta_{Z}^{i} \over \partial s_{\mu}}\right).}]$ The parameters are corrected by $[\delta s_{\mu}]$ and a new cycle of refinement is started until a minimum of E is reached. The weights $[w_{X} = 1/{\textstyle\sum\limits_{i=1}^{n}}(\Delta_{X}^{i})^{2}, \quad w_{Y} = 1/{\textstyle\sum\limits_{i=1}^{n}} (\Delta_{Y}^{i})^{2}, \quad w_{Z} = 1/{\textstyle\sum\limits_{i=1}^{n}} (\Delta_{Z}^{i})^{2}]$ are calculated with the current guess for $[s_{\mu}]$ at the beginning of each cycle.

The derivatives appearing in the normal equations can be worked out from the definitions given in Sections 11.3.2.2 and 11.3.2.4 , and only the form of the gradient of the Z residuals is shown. Assuming $[\sigma_{i} = \sigma_{M}/|\zeta_{i}|\ (i = 1,\ldots, n)]$ is constant for each reflection, the gradients of the Z residuals are obtained from the chain rule and the relation $[\hbox{d erf}(z)/\hbox{d}z = [2/(\pi)^{1/2}] \exp (-z^{2})]$ . $[\eqalign{{\partial \Delta_{Z}^{i} \over \partial s_{\mu}} &= {\partial \Delta_{Z}^{i} \over \partial \varphi_{i}} {\partial \varphi_{i} \over \partial s_{\mu}}\cr {\partial \Delta_{Z}^{i} \over \partial \varphi_{i}} &= {\Delta_{\varphi} \over (2\pi)^{1/2} \sigma_{i}} \sum\limits_{j=-\infty}^{\infty} \exp [-(\varphi_{0} + j\Delta_{\varphi} -\varphi_{i})^{2}/2\sigma_{i}^{2}]\cr {\partial \varphi_{i} \over \partial s_{\mu}} &= \cos \varphi_{i} {\partial \sin \varphi_{i} \over \partial s_{\mu}} - \sin \varphi_{i} {\partial \cos \varphi_{i} \over \partial s_{\mu}}.}]$ Obviously, $[\partial \Delta_{Z}^{i}/\partial s_{\mu}]$ is small for a fully recorded reflection because of the small values of all exponentials appearing in $[\partial \Delta_{Z}^{i}/\partial \varphi_{i}]$ . In contrast, the gradient for a partial reflection, equally recorded on two adjacent images, is most sensitive to parameter variations because one of the exponentials assumes its maximum value. In the limiting case of infinitely fine-sliced data, it can be shown that $[\lim_{\Delta_{\varphi \rightarrow 0}} \partial \Delta_{Z}^{i}/\partial \varphi_{i} = 1]$ . Thus, the refinement scheme based on observed Z centroids, as described here and implemented in XDS, is applicable to fine-sliced data – and to data recorded with a large oscillation range as well.

11.3.3. Integration

| top | pdf |

A fundamental requirement for a general integration method is that it should distinguish carefully between signal and background points within its integration domain. For weak reflections, this distinction cannot be made reliably because of the errors superimposed on the signal. The problem can be solved, however, provided that both weak and strong reflections share the same profile shape – an assumption that has been adopted by most data-processing packages.

The intensity distribution of a reflection can be modelled analytically or derived from the observed profiles of neighbouring strong spots. For the rotation method, the profile shape depends strongly on the specific path of the reflection through the Ewald sphere and on variations in the angle of incidence of the diffracted beam on a flat detector. These geometrical distortions can be eliminated by mapping the reflections onto the coordinate system defined in Section 11.3.2.3 , which simplifies the task of modelling the expected intensity distribution as all reflection profiles become similar.

11.3.3.1. Spot extraction

| top | pdf |

The region around a spot is defined by the two parameters $[\delta_{D}]$ and $[\delta_{M}]$ , which represent spot diameter and reflecting range, respectively. It is assumed that the coordinates of all image pixels contributing to the intensity of a spot satisfy $[|\varepsilon_{1}| \leq \delta_{D}/2, |\varepsilon_{2}| \leq \delta_{D}/2]$ and $[|\varepsilon_{3}| \leq \delta_{M}/2]$ when mapped to the profile coordinate system $[\{{\bf e}_{1}, {\bf e}_{2}, {\bf e}_{3}\}]$ defined in Section 11.3.2.3 . Regions of neighbouring reflections may overlap. As implemented in XDS, potential overlap is dealt with by a simple strategy: pixels within the overlap region are assigned to the nearest spot. This is carried out in two steps. First, reflections predicted to occur on a given rotation image are found by generating and testing all possible indices h, k, l up to the highest resolution recorded by the detector. Reflection indices, coordinates of the diffracted beam wave vector and the expected fraction of spot intensity recorded on the image are saved in a table. In the second step, each reflection boundary is traced in the image and corrected to exclude pixels belonging to overlapping reflections, which are rapidly located in the table by the hash technique. The image scaling factor obtained from the mean image background and the neighbourhood pixel values belonging to the reflections recorded in the image are saved on a scratch file dedicated to the currently processed data image.

At regular intervals, these files are merged such that all pixel values belonging to a spot found in the contributing images follow each other. Reflections for which contributing pixels are expected further ahead in data processing are just copied to a scratch output file. The other reflections are mapped to the Ewald sphere, as described below, and their three-dimensional profiles and accompanying information are routed to the main output file of the spot-extraction step. After the file-merging procedure, spot extraction continues.

11.3.3.2. Background

| top | pdf |

The region around a spot is assumed to have been chosen to be large enough to include a sufficient number of pixels which can be used for determination of the background. Background determination, as implemented in XDS, begins by sorting all pixels belonging to a reflection by increasing intensity. For weak or absent reflections, these values should represent a random sample drawn from a normal distribution. If this is not the case, the pixel with the largest intensity is removed until the sampling distribution of the remaining smaller items satisfies the expected distribution. This method will also exclude pixels with unexpectedly high values, such as ice reflections. The background, determined as the mean value of the accepted pixels, is systematically overestimated for strong spots because of some residual intensity extending into the accepted background pixels. This residual intensity is estimated from the expected distribution $[\omega (\varepsilon_{1}, \varepsilon_{2}, \varepsilon_{3})]$ defined in Section 11.3.2.3 and removed from the final background value.

11.3.3.3. Standard profiles

| top | pdf |

Reflection profiles are represented on the Ewald sphere within a domain $[D_{0}]$ comprising $[2n_{1} + 1, 2n_{2} + 1, 2n_{3} + 1]$ equidistant gridpoints along $[{\bf e}_{1}, {\bf e}_{2}, {\bf e}_{3}]$ , respectively. The sampling distances between adjacent grid points are then $[\Delta_{1} = \delta_{D}/(2n_{1} + 1)]$ , $[\Delta_{2} = \delta_{D}/(2n_{2} + 1), \Delta_{3} = \delta_{M}/(2n_{3} + 1)]$ . Thus, grid coordinate $[\nu_{3}\ (\nu_{3} = -n_{3},\ldots, n_{3})]$ covers the set of rotation angles $[\Gamma_{\nu_{3}} = \{\varphi' | (\nu_{3} - 1/2) \Delta_{3} \leq (\varphi' - \varphi)\cdot \zeta \leq (\nu_{3} + 1/2)\Delta_{3}\}.]$ Contributions to the spot intensity come from one or several adjacent data images $[(\;j = j_{1},\ldots,\ j_{2})]$ , each covering the set of rotation angles $[\Gamma_{j} = \{\varphi' | \varphi_{0} + (j - 1) \Delta_{\varphi} \leq \varphi' \leq \varphi_{0} + j\Delta_{\varphi}\}.]$ Assuming Gaussian profiles along $[{\bf e}_{3}]$ for all reflections (see Section 11.3.2.3 ), the fraction of counts (after subtraction of the background) contributed by data frame j to grid coordinate $[\nu_{3}]$ is $[\eqalign{f_{\nu_{3}j}&\approx {\textstyle\int\limits_{\Gamma_{j} \cap \Gamma_{\nu_{3}}}} \exp [-(\varphi' - \varphi)^{2}/2\sigma^{2}]\ \hbox{d} \varphi' \cr &\quad\times\bigg\{{\textstyle\int\limits_{\Gamma_{j}}} \exp [-(\varphi' - \varphi)^{2}/2\sigma^{2}]\ \hbox{d} \varphi'\bigg\}^{-1},}]$ where $[\sigma = \sigma_{M}/|\zeta|]$ . The integrals can be expressed in terms of the error function, for which efficient numerical approximations are available (Abramowitz & Stegun, 1972 ). Finally, each pixel on data image j belonging to the reflection is subdivided into $[5 \times 5]$ areas of equal size, and $[f_{\nu_{3}j}/25]$ of the pixel signal is added to the profile value at grid coordinates $[\nu_{1},\nu_{2},\nu_{3}]$ corresponding to each subdivision.

This complicated procedure leads to more uniform intensity profiles for all reflections than using their untransformed shape. This simplifies the task of modelling the expected intensity distribution needed for integration by profile fitting. As implemented in XDS, reference profiles are learnt every 5° of crystal rotation at nine positions on the detector, each covering an equal area of the detector face. In the learning phase, profile boxes of the strong reflections are normalized and added to their nearest reference profile boxes. The contributions are weighted according to the distance from the location of the reference profile. Each grid point within the average profile boxes is classified as signal if it is above 2% of the peak maximum. Finally, each profile is scaled such that the sum of its signal pixels normalizes to one. The analytic expression $[\omega(\varepsilon_{1},\varepsilon_{2},\varepsilon_{3})]$ defined in Section 11.3.2.3 for the expected intensity distribution is only a rough initial approximation which is now replaced by the empirical reference profiles.

11.3.3.4. Intensity estimation

| top | pdf |

If an expected intensity distribution $[\{p_{i} | i\in D_{0}\}]$ of the observed profile is given in a domain $[D_{0}]$ , the reflection intensity I can be estimated as $[I = {{\textstyle\sum\limits_{i\in D}}(c_{i}-b_{i})p_{i}/v_{i}}\bigg/{{\textstyle\sum\limits_{i\in D}}p_{i}^2/v_{i}},]$ which minimizes the function $[\psi(I) = {\textstyle\sum\limits_{i\in D}}(c_{i} - I\cdot p_{i} - b_{i})^2/v_{i} , \qquad {\textstyle\sum\limits_{i\in D_0}} p_{i} = 1.]$ $[b_{i},c_{i},v_{i}\ (i\in D)]$ are background, contents and variance of pixels observed in a subdomain $[D\subseteq D_{0}]$ of the expected distribution. The background $[b_{i}]$ underneath a diffraction spot is often assumed to be a constant which is estimated from the neighbourhood around the reflection. Determination of reflection intensities by profile fitting has a long tradition (Diamond, 1969 ; Ford, 1974 ; Kabsch, 1988b ; Otwinowski, 1993 ). Implementations of the method differ mainly in their assumptions about the variances $[v_{i}]$ . Ford uses constant variances, which works well for films, which have a high intrinsic background. In XDS, which was originally designed for a multiwire detector, $[v_{i}\propto p_{i}]$ was assumed, which results in a straight summation of background-subtracted counts within the expected profile region, $[I=\sum\nolimits_{i\in D}(c_{i}-b_{i})/\sum\nolimits_{i\in D}p_{i}]$ . This particular simple formula is very satisfactory for the low background typical of these detectors. For the general case, however, better results can be obtained by using $[v_{i} = b_{i} + Ip_{i}]$ for the pixel variances as shown by Otwinowski and implemented in DENZO and in the later version of XDS. Starting with $[v_{i} = b_{i}]$ , the intensity is now found by an iterative process which is terminated if the new intensity estimate becomes negative or does not change within a small tolerance, which is usually reached after three cycles. It can be shown that the solution thus obtained is unique.

11.3.4. Scaling

| top | pdf |

Usually, many statistically independent observations of symmetry-related reflections are recorded in the rotation images taken from one or several similar crystals of the same compound. The squared structure-factor amplitudes of equivalent reflections should be equal and the idea of scaling is to exploit this a priori knowledge to determine a correction factor for each observed intensity. These correction factors compensate to some extent for effects such as radiation damage, absorption, and variations in detector sensitivity and exposure times, as well as variations in size and disorder between different crystals.

The usual methods of scaling split the data into batches of roughly the same size, each covering one or more adjacent rotation images, and then determine a single scaling factor for all reflections in each batch. Neighbouring reflections may then receive quite different corrections if they are assigned to different batches. Since the selection of batch boundaries is to some extent arbitrary, a more continuous correction function would be preferable. This function could be modelled analytically (for example by using spherical harmonics) or empirically, as implemented in XSCALE and described below.

For each reflection, observational equations are defined as $[\psi_{hl\alpha} = (I_{hl} - g_{\alpha} I_{h})/\sigma_{hl}.]$ The subscript h represents the unique reflection indices and l enumerates all symmetry-related reflections to h. By definition, the unique reflection indices have the largest h, then k, then l value occuring in the set of all indices related by symmetry to the original indices, including Friedel mates. Thus, two reflections are symmetry-related if and only if their unique indices are identical. $[I_{h}]$ is the unknown `true' intensity and $[I_{hl}]$ , $[\sigma_{hl}]$ are symmetry-related observed intensities and their standard deviations, respectively. The subscript α denotes the coordinates at which the scaling function $[g_{\alpha}]$ should be evaluated. As implemented in XDS and XSCALE, $[\alpha = 1,\ldots,9]$ denotes nine positions uniformly distributed in the detector plane at the beginning of data collection, $[\alpha = 10,\ldots,18]$ the same positions on the detector but after the crystal has been rotated by, say, 5°, and so on. The scaling factors $[g_{\alpha}]$ and the estimated intensities $[I_{h}]$ are found at the minimum of the function $[\Psi = {\textstyle\sum\limits_{hl\alpha}}w_{hl\alpha}\Psi_{hl\alpha}^2.]$

The main difference from the method of Fox & Holmes (1966 ) is the introduction of the weights $[w_{hl\alpha}]$ . These weights depend upon the distance between each reflection hl and the positions α. They are monotonically decreasing functions of this distance, implemented as Gaussians in XDS and XSCALE. This results in a smoothing of the scaling factors since each reflection contributes to the observational equations in proportion to the weights $[w_{hl\alpha}]$ .

Minimization of Ψ is done iteratively. After each step, the $[g_{\alpha}]$ are replaced by $[\Delta g_{\alpha}+g_{\alpha}]$ and rescaled to a mean value of 1. The corrections $[\Delta g_{\alpha}]$ are determined from the normal equations $[{\textstyle\sum\limits_{\beta}}\ A_{\alpha\beta}\Delta g_{\beta}=b_{\alpha},]$ where $[\eqalign{A_{\alpha\beta} &= {\textstyle\sum\limits_h}[\delta_{\alpha\beta}I_h^2 u_{h\alpha} + (r_{h\alpha}v_{h\beta}+v_{h\alpha}r_{h\beta}-v_{h\alpha}v_{h\beta})/u_h]\cr b_{\alpha} &= {\textstyle\sum\limits_{h}} I_hr_{h\alpha}\cr I_{h} &= v_h/u_h\cr r_{h\alpha} &= v_{h\alpha}-g_{\alpha}u_{h\alpha}I_h\cr u_{h\alpha} &= {\textstyle\sum\limits_l} w_{hl\alpha}/\sigma_{hl}^2\cr v_{h\alpha} &= {\textstyle\sum\limits_l} w_{hl\alpha}I_{hl}/\sigma_{hl}^2\cr u_h &= {\textstyle\sum\limits_\alpha} g_{\alpha}^2u_{h\alpha}\cr v_h &= {\textstyle\sum\limits_\alpha} g_{\alpha}v_{h\alpha}.}]$

In case a `true' intensity $[I_{h}]$ is available from a reference data set, the non-diagonal elements are omitted from the sum over h in the normal matrix $[A_{\alpha\beta}]$ . The corrections $[\Delta g_{\alpha}]$ are expanded in terms of the eigenvectors of the normal matrix, thereby avoiding shifts along eigenvectors with very small eigenvalues (Diamond, 1966 ). This filtering method is essential since the normal matrix has zero determinant if no reference data set is available.

11.3.5. Post refinement

| top | pdf |

The number of fully recorded reflections on each single image rapidly declines for small oscillation ranges and the complete intensities of the partially recorded reflections have to be estimated. This presented a serious obstacle in early structural work on virus crystals, as the crystal had to be replaced after each exposure on account of radiation damage. A solution of this problem, the `post refinement' technique, was found by Schutt, Winkler and Harrison, and variants of this powerful method have been incorporated into most data-reduction programs [for a detailed discussion, see Harrison et al. (1985 ); Rossmann (1985 )]. The method derives complete intensities of reflections only partially recorded on an image from accurate estimates for the fractions of observed intensity, the `partiality'. The partiality of each reflection can always be calculated as a function of orientation, unit-cell metric, mosaic spread of the crystal and model intensity distributions. Obviously, the accuracy of the estimated full reflection intensity then strongly depends on a precise knowledge of the parameters describing the diffraction experiment. Usually, for many of the partial reflections, symmetry-related fully recorded ones can be found, and the list of such pairs of intensity observations can be used to refine the required parameters by a least-squares procedure. Clearly, this refinement is carried out after all images have been processed, which explains why the procedure is called `post refinement'.

Adjustments of the diffraction parameters $[s_{\mu}\ (\mu = 1,\ldots,k)]$ are determined by minimization of the function E, which is defined as the weighted sum of squared residuals between calculated and observed partial intensites. $[\eqalign{E &= {\textstyle\sum\limits_{hj}} w_{hj} (\Delta_{hj})^2\cr \Delta_{hj} &= R_{j} (\varphi_{hj})g_{j} I_{h} - I_{hj}\cr w_{hj} &= 1/ \{\sigma^{2}(I_{hj}) + [R_{j}(\varphi_{hj})g_{j}]^{2} \sigma^{2}(I_{h})\}.}]$ Here, $[I_{hj}]$ is the intensity recorded on image j of a partial reflection with indices summarized as $[hj,\ I_{h}]$ is the mean of the observed intensities of all fully recorded reflections symmetry-equivalent to $[hj,\ g_{j}]$ is the inverse scaling factor of image $[j, \varphi_{hj}]$ is the calculated spindle angle of reflection hj at diffraction and $[R_{j}]$ is the computed fraction of total intensity recorded on image j.

Expansion of the residuals $[\Delta_{hj}]$ to first order in the parameter changes $[\delta s_{\mu}]$ and minimization of $[E(\delta s_{\mu})]$ leads to the k normal equations $[\sum\limits_{\mu'=1}^{k} \left(\sum\limits_{hj} w_{hj} {\partial\Delta_{hj} \over \partial s_\mu} {\partial\Delta_{hj} \over \partial s_{\mu'}}\right) \delta s_{\mu'} = - \sum\limits_{hj} w_{hj}\Delta_{hj} {\partial\Delta_{hj} \over \partial s_{\mu}}.]$ Often, the normal matrix is ill-conditioned, since changes in some unit-cell parameters or small rotations of the crystal about the incident X-ray beam do not significantly affect the calculated partiality $[R_{j}]$ . To take care of these difficulties, the system of equations is rescaled to yield unit diagonal elements for the normal matrix and the correction vector $[\delta s_{\mu}]$ is filtered by projection into a subspace defined by the eigenvectors of the normal matrix with sufficiently large eigenvalues (Diamond, 1966 ).

The parameters are corrected by the filtered $[\delta s_{\mu}]$ and a new cycle of refinement is started until a minimum of E is reached. The weights, residuals and their gradients are calculated using the current values for $[s_{\mu}]$ and $[g_{j}]$ at the beginning of each cycle. The derivatives $[{\partial\Delta_{hj} \over \partial s_{\mu}} = g_{j} I_{h} \left({\partial R_{j} \over \partial\varphi_{hj}} {\partial \varphi_{hj} \over \partial s_{\mu}} + {\partial R_{j} \over \partial\sigma_{M}} {\partial\sigma_{M} \over \partial s_{\mu}} + {\partial R_{j} \over \partial|\zeta_{hj}|} {\partial|\zeta_{hj}| \over \partial s_{\mu}}\right)]$ appearing in the normal equations can be worked out from the definitions given in Sections 11.3.2.2 and 11.3.2.4 (to simplify the following equations, the subscript hj is omitted). The fraction $[R_{j}]$ of the total intensity can be expressed in terms of the error function (see Section 11.3.2.4 ) as $[\eqalign{R_{j} &= [\hbox{erf} (z_{1}) - \hbox{erf}(z_{2})]/2\cr z_{1} &= |\zeta|(\varphi_{0} + j\Delta_{\varphi} - \varphi)/(2)^{1/2}\sigma_{M}\cr z_{2} &= |\zeta| [\varphi_{0} + (\;j - 1)\Delta_{\varphi} - \varphi]/(2)^{1/2}\sigma_{M}.}]$ Using the relation $[\hbox{d erf}(z)/\hbox{d}z = [2/(\pi)^{1/2}] \exp(-z^{2})]$ , the derivatives of $[R_{j}]$ are $[\eqalign{\partial R_{j}/\partial\varphi &= [\exp(- z_{2}^{2}) - \exp (- z_{1}^{2})] |\zeta|/[\sigma_{M} (2\pi)^{1/2}]\cr \partial R_{j}/\partial\sigma_{M} &= [z_{2} \exp (- z_{2}^{2}) - z_{1} \exp(- z_{1}^{2})]/[\sigma_{M} (\pi)^{1/2}]\cr \partial R_{j}/\partial|\zeta| &= [z_{1} \exp (- z_{1}^{2}) - z_{2} \exp (- z_{2}^{2})]/[|\zeta| (\pi)^{1/2}].}]$ It remains to work out the derivatives $[\partial \varphi/\partial s_{\mu}]$ , $[\partial\sigma_{M}/\partial s_{\mu}]$ and $[\partial |\zeta|/\partial s_{\mu}]$ (not shown here). As discussed in detail by Greenhough & Helliwell (1982 ), spectral dispersion and asymmetric beam cross fire lead to some variation of $[\sigma_{M}]$ , which makes it necessary to include additional parameters in the list $[s_{\mu}]$ . The effect of these parameters on the partiality is dealt with easily by the derivatives $[\partial\sigma_{M}/\partial s_{\mu}]$ .

The refinement scheme described above requires initial scaling factors $[g_{j}]$ . With the now improved estimates for the partialities $[R_{j}]$ , a new set of scaling factors can be obtained by the method outlined in Section 11.3.4 . This alternating procedure of scaling and post refinement usually converges within three cycles.

The use of error functions for modelling partiality, as implicated by a Gaussian model for describing spot shape, was chosen here for reasons of conceptual simplicity and coherence. This choice is unlikely to alter significantly the results of post refinement that are based on other functions of similar form [see the discussion by Rossmann (1985 )].

11.3.6. Space-group assignment

| top | pdf |

Identification of the correct space group is not always an easy task and should be postponed for as long as possible. Fortunately, all data processing as implemented in the program XDS can be carried out even in the absence of any knowledge of crystal symmetry and cell constants. In this case, a reduced cell is extracted from the observed diffraction pattern and processing of the data images continues to completion as if the crystal were triclinic. Clearly, the reflection indices then refer to the reduced cell and must be reindexed once the space group is known. For all space groups, the required reindexing transformation is linear and involves only whole numbers as shown in Part 9 of IT A. The following description and example are taken from Kabsch (1993 ).

Space-group assignment is carried out in two steps under control of the crystallographer once integrated intensities of all reflections are available. First, the Bravais lattices that are compatible with the observed reduced cell are identified. In the second step, any of the plausible space groups may be tested and rated according to symmetry R factors and systematic absences of integrated reflection intensities after reindexing. Additional acceptance criteria are obtained from refinement, now using a reduced set of independent parameters describing the conventional unit cell which should not lead to a significant increase of r.m.s. deviations between observed and calculated reflection positions and angles.

11.3.6.1. Determination of the Bravais lattice

| top | pdf |

The determination of possible Bravais lattices is based upon the concept of the reduced cell whose metric parameters characterize 44 lattice types as described in Part 9 of IT A. A primitive basis $[{\bf b}_{1}, {\bf b}_{2}, {\bf b}_{3}]$ of a given lattice is defined there as a reduced cell if it is right-handed and if the components of its metric tensor $[\matrix{A = {\bf b}_{1}\cdot {\bf b}_{1}, &\qquad B = {\bf b}_{2}\cdot {\bf b}_{2}, &\qquad C = {\bf b}_{3}\cdot {\bf b}_{3},\cr D = {\bf b}_{2}\cdot {\bf b}_{3}, &\qquad E = {\bf b}_{1}\cdot {\bf b}_{3}, &\qquad F = {\bf b}_{1}\cdot {\bf b}_{2}\cr}]$ satisfy a number of conditions (inequalities). The main conditions state that the basis vectors are the shortest three linear independent lattice vectors with either all acute or all non-acute angles between them. As specified in IT A, each of the 44 lattice types is characterized by additional equality relations among the six components of the reduced-cell metric tensor. As an example, for lattice character 13 (Bravais type oC) the components of the metric tensor of the reduced cell must satisfy $[A = B, \quad B \leq C, \quad D = 0, \quad E = 0, \quad 0 \leq - F \leq A/2.]$

Any primitive triclinic cell describing a given lattice can be converted into a reduced cell. It is well known, however, that the reduced cell thus derived is sensitive to experimental error. Hence, the direct approach of first deriving the correct reduced cell and then finding the lattice type is unstable and may in certain cases even prevent the identification of the correct Bravais lattice.

A suitable solution of the problem has been found that avoids any decision about what the `true' reduced cell is. The essential requirements of this procedure are: (a) a database of possible reduced cells and (b) a backward search strategy that finds the best-fitting cell in the database for each lattice type.

The database is derived from a seed cell which strictly satisfies the definitions for a reduced cell. All cells of the same volume as the seed cell whose basis vectors can be linearly expressed in terms of the seed vectors by indices $[-1,\ 0,]$ or +1 are included in the database. Each unit cell in the database is considered as a potential reduced cell even though some of the defining conditions as given in Part 9 of IT A may be violated. These violations are treated as being due to experimental error.

The backward search strategy starts with the hypothesis that the lattice type is already known and identifies the best-fitting cell in the database of possible reduced cells. Contrary to a forward directed search, it is now always possible to decide which conditions have to be satisfied by the components of the metric tensor of the reduced cell. The total amount by which all these equality and inequality conditions are violated is used as a quality index. This measure is defined below for lattice type 13 oC testing a potential reduced cell $[{\bf b}_{1}, {\bf b}_{2}, {\bf b}_{3}]$ from the database for agreement. Positive values of the quality index $[p_{13}]$ indicate that some conditions are not satisfied. $[\eqalign{p_{13}({\bf b}_{1}, {\bf b}_{2}, {\bf b}_{3}) &= |A - B| + \max(0,B - C) + |D| + |E|\cr &\quad+ \max(0,F) + \max(0, - F - A/2).}]$ All potential reduced cells in the database are tested and the smallest value for $[p_{13}]$ is assigned to lattice type 13. This test is carried out for all 44 possible lattice types using quality indices derived in a similar way from the defining conditions as listed in Part 9 of IT A. For each of the 44 lattice types thus tested, the procedure described here returns the quality index, the conventional cell parameters and a transformation matrix relating original indices with respect to the seed cell to the new indices with respect to the conventional cell. These index-transformation matrices are derived from those given in Table 9.2.5.1 in IT A.

The results obtained by this method are shown in Table 11.3.6.1 for the example of a 1.5° oscillation data film containing 1313 strong diffraction spots which were located automatically. The space group of the crystal is $[C222_{1}]$ and the cell constants are $[a = 72.9, b = 100.1, c = 92.6\;\hbox{\AA}]$ . The entry for the correct Bravais lattice oC with derived cell constants close to the true ones has a low value for its quality index and thus appears as a possible explanation of the observed diffraction pattern.

Table 11.3.6.1| top | pdf |
Rating of lattice types implied by a given reduced cell

Lattice type	Quality index	Conventional cell constants (Å, °)						Reindexing transformation
Lattice type	Quality index	a	b	c	α	β	γ	Reindexing transformation
1 cF	999.0	119.3	137.3	119.1	121.1	77.5	122.6	$[11\bar{1}0/1\bar{1}10/\bar{1}\bar{1}\bar{1}0]$
2 hR	770.1	74.6	111.7	137.4	103.6	89.1	108.8	$[1100/\bar{1}0\bar{1}0/\bar{1}110]$
3 cP	769.7	62.1	63.5	92.9	90.0	90.1	107.2
5 cI	936.0	111.7	74.6	112.6	70.1	53.6	71.2
4 hR	769.5	101.1	111.9	119.1	116.5	89.1	116.4	$[1\bar{1}00/\bar{1}010/\bar{1}\bar{1}\bar{1}0]$
6 tI	999.0	112.6	111.7	74.6	71.2	70.1	53.6
7 tI	999.0	111.7	74.6	112.6	70.1	53.6	71.2
8 oI	999.0	74.6	111.7	112.6	53.6	70.1	71.2	$[\bar{1}\bar{1}00/\bar{1}0\bar{1}0/0\bar{1}\bar{1}0]$
9 hR	772.7	62.1	74.6	296.6	90.5	105.8	125.5	$[1000/\bar{1}\bar{1}00/\bar{1}1\bar{3}0]$
10 mC	24.0	101.1	74.6	92.9	90.1	90.0	91.3	$[1\bar{1}00/1100/0010]$
11 tP	174.8	62.1	63.5	92.9	90.0	90.1	107.2
12 hP	122.8	62.1	63.5	92.9	90.0	90.1	107.2
13 oC	23.8	74.6	101.1	92.9	90.0	90.1	88.7	$[1100/\bar{1}100/0010]$
15 tI	672.7	62.1	63.5	200.2	77.0	77.6	107.2
16 oF	999.0	74.6	101.1	200.2	90.5	111.8	88.7	$[\bar{1}\bar{1}00/1\bar{1}00/1120]$
14 mC	23.4	74.6	101.1	92.9	90.0	90.1	88.7	$[1100/\bar{1}100/0010]$
17 mC	999.0	101.1	74.6	111.7	71.2	116.4	88.7	$[1\bar{1}00/\bar{1}\bar{1}00/\bar{1}0\bar{1}0]$
18 tI	999.0	112.6	119.1	62.1	68.7	99.5	115.4	$[01\bar{1}0/1110/1000]$
19 oI	999.0	62.1	112.6	119.1	64.6	68.7	80.5	$[\bar{1}000/01\bar{1}0/\bar{1}\bar{1}\bar{1}0]$
20 mC	746.3	112.6	112.6	62.1	99.5	99.6	111.3	$[0\bar{1}\bar{1}0/0\bar{1}10/\bar{1}000]$
21 tP	748.0	63.5	92.9	62.1	90.1	107.2	90.0
22 hP	999.0	63.5	92.9	62.1	90.1	107.2	90.0
23 oC	747.8	112.6	112.6	62.1	80.5	99.6	68.7	$[0110/0\bar{1}10/1000]$
24 hR	999.0	154.8	112.6	62.1	80.5	80.9	84.3	$[1210/0\bar{1}10/1000]$
25 mC	746.1	112.6	112.6	62.1	80.5	99.6	68.7	$[0110/0\bar{1}10/1000]$
26 oF	624.9	62.1	123.9	195.9	86.4	108.4	101.5	$[1000/\bar{1}\bar{2}00/\bar{1}0\bar{2}0]$
27 mC	499.7	123.9	62.1	112.6	80.5	119.7	78.5	$[\bar{1}\bar{2}00/\bar{1}000/01\bar{1}0]$
28 mC	325.0	62.1	195.9	63.5	95.4	107.2	71.6	$[\bar{1}000/\bar{1}0\bar{2}0/0\bar{1}00]$
29 mC	99.8	62.1	123.9	92.9	90.0	90.1	78.5
30 mC	336.4	63.5	196.5	62.1	95.4	107.2	71.1	$[0\bar{1}00/0\bar{1}20/\bar{1}000]$
31 aP	0.2	62.1	63.5	92.9	90.0	89.9	72.8	$[1000/0\bar{1}00/00\bar{1}0]$
32 oP	152.0	62.1	63.5	92.9	90.0	90.1	107.2
40 oC	413.0	63.5	196.4	62.1	84.5	107.2	108.9	$[0\bar{1}00/0120/\bar{1}000]$
35 mP	151.8	63.5	62.1	92.9	90.1	90.0	107.2	$[0\bar{1}00/\bar{1}000/00\bar{1}0]$
36 oC	400.3	62.1	195.9	63.5	84.6	107.2	108.4	$[1000/\bar{1}0\bar{2}0/0100]$
33 mP	151.2	62.1	63.5	92.9	90.0	90.1	107.2
38 oC	100.1	62.1	123.9	92.9	90.0	90.1	101.5	$[\bar{1}000/1200/00\bar{1}0]$
34 mP	1.0	62.1	92.9	63.5	90.0	107.2	90.1	$[\bar{1}000/00\bar{1}0/0\bar{1}00]$
42 oI	661.3	62.1	63.5	200.2	103.0	102.4	107.2	$[\bar{1}000/0\bar{1}00/1120]$
41 mC	412.2	196.4	63.5	62.1	107.2	95.5	71.1	$[0\bar{1}\bar{2}0/0\bar{1}00/\bar{1}000]$
37 mC	400.1	195.9	62.1	63.5	107.2	95.4	71.6
39 mC	99.9	123.9	62.1	92.9	90.1	90.0	78.5	$[\bar{1}\bar{2}00/\bar{1}000/00\bar{1}0]$
43 mI	999.0	74.6	200.2	63.5	103.0	127.3	68.2	$[1100/1120/0\bar{1}00]$
44 aP	0.0	62.1	63.5	92.9	90.0	90.1	107.2

11.3.6.2. Finding possible space groups

| top | pdf |

Inspection of the table rating the likelihood of each of the 44 lattice types usually reveals a rather limited set of possible space groups. Furthermore, the absence of parity-changing symmetry operators required for protein crystals restricts the number of possible space groups to 65 instead of 230. Any space group can be tested by repeating only the final steps of data processing. These steps include a comparison of symmetry-related reflection intensities, as well as a refinement of the parameters controlling the diffraction pattern after reindexing the reflections by the appropriate transformation. Low r.m.s. deviations between the observed and refined spot positions, as well as small R factors for symmetry-related reflection intensities, indicate that the constraints imposed by the tentatively chosen space group are satisfied. The space group with highest symmetry compatible with the data is almost certainly correct if the data set is sufficiently complete and redundant, which requires that each symmetry element relates a sufficient number of reflections to one another.

For the example of a 1.5° oscillation data film given above, space-group determination consists of the following steps. Inspection of Table 11.3.6.1 indicates that lattice characters 10, 13, 14 and 34, besides the triclinic characters 31 and 44, are approximately compatible with the observed diffraction pattern. The highest lattice symmetry is orthorhombic (character 13, Bravais type oC), which limits the possible space groups for protein crystals to either $[C{222}_{1}]$ or C222. Processing of all films in the data set was completed in space group P1 using the cell constants shown for lattice character 44. To test whether the crystal has space-group symmetry C222 and conventional cell constants $[a = 74.6, b = 101.1, c = 92.9\;\hbox{\AA}]$ , the final steps of data processing were repeated after reindexing the reflections by the transformation $[h' = 1\cdot h + 1\cdot k + 0\cdot l + 0]$ , $[k' = -1\cdot h + 1\cdot k + 0\cdot l + 0]$ , $[l' = 0\cdot h + 0\cdot k + 1\cdot l + 0]$ as specified for lattice character 13. Note that the transformation also provides a simple tool for correcting the indices if all reflections are misindexed by a constant. The results clearly show that the crystal has space-group symmetry $[C222_{1}]$ . The presence of the $[2_{1}]$ axis was deduced from the rather weak intensities observed for reflections of type $[0 0 l' = \hbox{odd}]$ .

References

Abramowitz, M. & Stegun, I. A. (1972). Handbook of mathematical functions. New York: Dover Publications.Google Scholar

Bricogne, G. (1986). Indexing and the Fourier transform. In Proceedings of the EEC cooperative workshop on position-sensitive detector software (phase III), p. 28. LURE, 12–19 November.Google Scholar

Diamond, R. (1966). A mathematical model-building procedure for proteins. Acta Cryst. 21, 253–266.Google Scholar

Diamond, R. (1969). Profile analysis in single crystal diffractometry. Acta Cryst. A25, 43–55.Google Scholar

Dijkstra, E. W. (1976). A discipline of programming, pp. 154–167. New Jersey: Prentice-Hall.Google Scholar

Ford, G. C. (1974). Intensity determination by profile fitting applied to precession photographs. J. Appl. Cryst. 7, 555–564.Google Scholar

Fox, G. C. & Holmes, K. C. (1966). An alternative method of solving the layer scaling equations of Hamilton, Rollett and Sparks. Acta Cryst. 20, 886–891.Google Scholar

Greenhough, T. J. & Helliwell, J. R. (1982). Oscillation camera data processing: reflecting range and prediction of partiality. I. Conventional X-ray sources. J. Appl. Cryst. 15, 338–351.Google Scholar

Harrison, S. C., Winkler, F. K., Schutt, C. E. & Durbin, R. M. (1985). Oscillation method with large unit cells. Methods Enzymol. 114A, 211–237.Google Scholar

Howard, A. (1986). Autoindexing. In Proceedings of the EEC cooperative workshop on position-sensitive detector software (phases I & II), pp. 89–94. LURE, 26 May–7 June.Google Scholar

International Tables for Crystallography (2005). Vol. A. Space-group symmetry, edited by Th. Hahn, ch. 9.2. Heidelberg: Springer.Google Scholar

Kabsch, W. (1988a). Automatic indexing of rotation diffraction patterns. J. Appl. Cryst. 21, 67–71.Google Scholar

Kabsch, W. (1988b). Evaluation of single-crystal X-ray diffraction data from a position-sensitive detector. J. Appl. Cryst. 21, 916–924.Google Scholar

Kabsch, W. (1993). Automatic processing of rotation diffraction data from crystals of initially unknown symmetry and cell constants. J. Appl. Cryst. 26, 795–800.Google Scholar

Otwinowski, Z. (1993). Oscillation data reduction program. In Proceedings of the CCP4 study weekend. Data collection and processing, edited by L. Sawyer, N. Isaacs & S. Bailey, pp. 56–62. Warrington: Daresbury Laboratory.Google Scholar

Otwinowski, Z. & Minor, W. (1997). Processing of X-ray diffraction data collected in oscillation mode. Methods Enzymol. 276, 307–326.Google Scholar

Pflugrath, J. W. (1997). Diffraction-data processing for electronic detectors: theory and practice. Methods Enzymol. 276A, 286–306.Google Scholar

Rossmann, M. G. (1985). Determining the intensity of Bragg reflections from oscillation photographs. Methods Enzymol. 114A, 237–280.Google Scholar

Schutt, C. & Winkler, F. K. (1977). The oscillation method for very large unit cells. In The rotation method in crystallography, edited by U. W. Arndt & A. J. Wonacott, pp. 173–186. Amsterdam, New York, Oxford: North-Holland.Google Scholar

Steller, I., Bolotovsky, R. & Rossmann, M. G. (1997). An algorithm for automatic indexing of oscillation images using Fourier analysis. J. Appl. Cryst. 30, 1036–1040.Google Scholar

Wirth, N. (1976). Algorithms + data structures = programs, pp. 264–274. New York: Prentice-Hall.Google Scholar

International Tables for Crystallography (2006). Vol. F. ch. 11.3, pp. 218-225
https://doi.org/10.1107/97809553602060000676