A 3D profile function suitable for integration of neutron time-of-flight single crystal diffraction peaks

A 3D profile function suitable for integration of neutron time-of-flight single crystal diffraction peaks

Nuclear Instruments and Methods in Physics Research A (xxxx) xxxx–xxxx Contents lists available at ScienceDirect Nuclear Instruments and Methods in ...

1MB Sizes 3 Downloads 156 Views

Nuclear Instruments and Methods in Physics Research A (xxxx) xxxx–xxxx

Contents lists available at ScienceDirect

Nuclear Instruments and Methods in Physics Research A journal homepage: www.elsevier.com/locate/nima

Technical notes

A 3D profile function suitable for integration of neutron time-of-flight single crystal diffraction peaks Matthias J. Gutmann Rutherford Appleton Laboratory, ISIS Facility, Chilton Didcot, Oxfordshire OX11 0QX, United Kingdom

A BS T RAC T A 3D profile function is presented suitable to integrate reflections arising in time-of-flight (TOF) single crystal neutron diffraction experiments. In order to account for the large asymmetry of the peak shape in the TOF direction, a 3D Gaussian ellipsoid in the pixel (x, z) and time-of-flight coordinates is convoluted with a rising and falling exponential along the time-of-flight direction. An analytic expression is derived, making it suitable for least-squares fitting. The application of this function in detector space or reciprocal space is straightforward.

1. Introduction The integration of Bragg peaks is among the key steps in the data reduction of single crystal diffraction experiments. At a spallation neutron source, when pixelated area detectors are used, each pixel has time-of-flight resolution, i.e. it records the arrival time of every neutron after emission from the spallation target/moderator assembly. The interaction of the neutrons within the moderator leads to an asymmetric peak shape along the time-of-flight direction (Fig. 1). Examples of this type of instruments are SXD at ISIS (Oxfordshire, United Kingdom), TOPAZ, CORELLI and MANDI at SNS (Oak Ridge National Laboratory, USA) or SENJU and IBIX at J-PARC (Ibaraki prefecture, Japan) [1–6]. Hence, single crystal reflections are not perfect 3D Gaussian ellipsoids but typically have a steep rising edge on the low TOF side and an elongated tail on the high TOF side. This is also the case in TOF neutron powder experiments. Various approaches to integrating such peaks have been published in the literature and shown to give good results. A peak shape independent approach is the seed-skewness method originally developed for the case of 2D X-ray images and extended to 3D for the TOF case [7,8]. Profile fitting methods based on 1D fitting are widely used, in which the peak is fitted with a 1D profile in each pixel and the results summed. Various profile methods have been proposed to capture the asymmetric tails, e.g. Vavilov, Landau convoluted with a Gaussian and GSAS type I functions [9–12]. The latter appears the most popular and is also extensively used in neutron TOF powder diffraction. This method can be applied in either detector space or reciprocal space. In detector space, the peaks are integrated with respect to the detector coordinates which are the (x, z) pixel and TOF coordinates. For reciprocal space integration, the (x, z, t) coordinates are converted to reciprocal space (Qx, Qy, Qz). In this case, the asymmetry along TOF is translated into an asymmetry which is pointing radially from the peak position to the origin of reciprocal

space. To our knowledge, no 3D profile function suitable to integrate single crystal peaks from TOF neutron instruments has been published to date. In this paper an account is given of the 3D profile fitting function that has been implemented in the SXD2001 package and has been in use since 2004 on the SXD instrument [13]. This function not only gives integrated intensities for crystal structure solution/refinement but also provides peak positions that can be used to reliably refine lattice parameters. Although the emphasis is on using this function with a 2D pixelated detector due to the context in which it was developed, the function can generally be applied where such a peak shape arises. In the following the profile function is presented in detail together with an extension to its use in reciprocal space directly. A speed comparison between a serial and parallel version of the computational routine and some outlook on further uses of this function are also given. 2. Profile function In time-of-flight single-crystal neutron diffraction using a timeresolved, pixelated area detector one can describe the peak shape of Bragg reflections by a three-dimensional Gaussian ellipsoid along the (x,z) pixel and time-of-flight coordinates. The asymmetry along timeof-flight is taken into account by convoluting this function with a backto-back exponential. In the following, the peak shape shall be derived to a closed expression. The starting point is a three-dimensional Gaussian ellipsoid described as follows: →T →

G (Δx, Δz , Δt ) = e−1/2 v H v (1) → Here, v is a column vector and H a symmetric 3×3 matrix given by:

http://dx.doi.org/10.1016/j.nima.2016.12.026 Received 29 September 2016; Received in revised form 9 November 2016; Accepted 14 December 2016 0168-9002/ © 2016 Elsevier B.V. All rights reserved.

Please cite this article as: Gutmann, M.J., Nuclear Instruments and Methods in Physics Research A (2016), http://dx.doi.org/10.1016/j.nima.2016.12.026

Nuclear Instruments and Methods in Physics Research A (xxxx) xxxx–xxxx

M.J. Gutmann


α + H6 Δt + H3 Δx + H5 Δz 2H6


β − H6 Δt − H3 Δx − H5 Δz 2H6

Erfc (x ) is the complementary Gaussian error function given by:

Erfc (x ) =

2 π



e−t dt

Additional information about the resolution along TOF as a function of scattering angle 2θ can be incorporated and proves advantageous as a starting value for peak fitting and in the integration of weak reflections:

⎛ dt ⎞2 ⎛ dL ⎞2 dt = (dθ ∙cotan (θ ))2 + ⎜ 0 ⎟ + ⎜ ⎟ ⎝ L ⎠ ⎝ t0 ⎠ t

Here dθ is the angle dependent constant of the diffractometer resolution, dL/L the uncertainty in flight path and dt0/t0 the intrinsic pulse width from the moderator which can be obtained from a peak close to backscattering or the known moderator characteristics. Eq. (7) multiplied with the TOF, inverted and squared would correspond to the matrix element H6. The normalisation of Eq. (6) can be carried out numerically. Alternatively, the convolution along the time-of-flight direction can be carried out using a fast Fourier transform, in which case it becomes a multiplication in Fourier space. Fig. 2 illustrates a fit to a TOF single crystal reflection from a NaCl sample using this 3D profile function and a linear background in x, z, and t. In particular it should be noted that the peak shape asymmetry is captured quite well. Although the function presented here might seem computationally expensive, particularly in the use of exponentials, (profile fitting can encompass several thousand points simultaneously), its analytic form is well suited for parallelisation with respect to the evaluation of a given set of discrete data points. Test runs of a serial and parallel version in C++ and OpenCL over 10,000 function evaluations indicate an over 70 fold speed increase of the parallel version run on 32 cores (~0.09 s) compared to a single core of the same CPU (~7.3 s). To set this into context, fitting the peak shown in Fig. 2 required 657 function calls. This would mean that 10,000 peaks could be fitted in about 1 min on 32 CPU cores.

Fig. 1. A typical peak profile recorded on a pixel in the backscattering detector module of the SXD instrument from a single crystal of NaCl corresponding to the (2, −6,0) reflection. The line connecting the points is a guide to the eye.

⎛ H1 H2 H3 ⎞ ⎛ Δx ⎞ ⎛ x − x 0 ⎞ → v = ⎜⎜ Δz ⎟⎟ = ⎜⎜ z − z 0 ⎟⎟ and H = ⎜⎜ H2 H4 H5 ⎟⎟ ⎝ Δt ⎠ ⎝ t − t0 ⎠ ⎝ H3 H5 H6 ⎠ x 0 , z 0 , and t0 describe the peak centroid with respect to the pixel and time-of-flight coordinates. H describes the “inverse squared width”. The off-diagonal elements allow for an arbitrary orientation of the ellipsoid on the detector. Eq. (1) may be expanded to yield: G (Δx, Δz, Δt ) = e−1/2(H1Δx

2 +2H Δx Δz +2H Δx Δt + H Δz 2 +2H ΔzΔt + H Δt 2 ) 2 3 4 5 6


The normalisation factor for such a function has been derived by Sivia in [14] and is given by:

det(H ) (2π )3/2

NG =


ensuring that the Gaussian ellipsoid has a volume of 1.0. The asymmetry mainly shows up along the time-of-flight direction. This can be accurately modelled by convoluting the Gaussian ellipsoid (1) with a back-to-back exponential. The resulting peak shape is then:

⎧ 0 ⎪ ∫ G (Δx, Δz, Δt − τ )⋅2N ⋅e ατ dτ , τ < 0 H (Δx, Δz , Δt ) = ⎨ −∞ ∞ ⎪ ∫ G (Δx, Δz, Δt − τ )⋅2N ⋅e−βτ dτ , τ > 0 ⎩ 0

3. Summary and conclusion In this paper, a 3D profile function is presented as a natural extension of a 1D peak profile based on a Gaussian ellipsoid convoluted with a rising and falling exponential accounting for the asymmetry along TOF. The 1D Gaussian is replaced with a 3D Gauss ellipsoid. A formula is obtained in closed form making it suitable for least squares refinement. The integration not only returns integrated intensities for solving and refining crystal structures but also peak positions leading to lattice parameters which are consistent with those from neutron TOF powder diffraction and laboratory X-rays [15]. There are also advantages in integrating weak reflections. This often means that parameters of a function need to be constrained or fewer parameters refined. The advantage here is that such constraints can be used on the width in all three dimensions, i.e. the H matrix, α and β. These can be constrained from nearby strong reflections and using Eq. (7). This is not so straightforward to do for 1D or 2D based integrations. A similar approach applies for partially overlapping peaks. The returned ellipsoid parameters can also be used to compare widths of reflections of various samples, e.g. for strain comparisons, similar as was done in [16]. The function can be readily applied in Q space where (x, z, t) is replaced by (Qx, Qy, Qz) and likewise for the inverse squared variance matrix H [17]. A parallel version of the routine is shown to lead to a dramatic speed-up of the function evaluation. This opens up the possibility of using this function directly as an option for peak searching which means integrated intensities could be obtained as one of the very first steps in the data processing rather than one of the last steps.


α and β are the inverse time constants of the rising and falling exponential, respectively. The normalisation factor of the convolution is given by: 2N =

αβ (α + β )


After some laborious calculations the resulting peak shape function can be expressed in a closed form as follows:

H (Δx, Δz, Δt ) = 2N ⋅NG⋅

π − 1 H Δx 2 − H2 Δx Δz − 1 H4 Δz 2 + 2 ⋅e 2 1 2H6

⋅(eu⋅Erfc ( y) + e v⋅Erfc (w ))

H32 Δx 2 +2H3 H5 Δx Δz + H52 Δz 2 2H6




1 α (α + 2H6 Δt + 2H3 Δx + 2H5 Δz ) 2


1 β (β − 2H6 Δt − 2H3 Δx − 2H5 Δz ) 2



Nuclear Instruments and Methods in Physics Research A (xxxx) xxxx–xxxx

M.J. Gutmann

Fig. 2. Fitted profile of the reflection shown in Fig. 1 from a single crystal of NaCl. The data are shown in a checkerboard pattern corresponding to an array of 7×7 (x, z) pixels containing the peak. Each checkerboard field contains the TOF range around the peak that has been fitted. They grey line is the fit. Note that the intensity scale is separate for each pixel and is to the left in each checkerboard field.

[3] L. Coates, M.J. Cuneo, M.J. Frost, J. He, K.L. Weiss, S.J. Tomanicek, H. McFeeters, V.G. Vandavasi, P. Langan, E.B. Iverson, J. Appl. Crystallogr. 48 (2015) 1302. [4] S. Rosenkranz, R. Osborn, Pramana- J. Phys. 71 (2008) 705. [5] T. Ohhara, R. Kiyanagi, K. Oikawa, K. Kaneko, T. Kawasaki, I. Tamura, A. Nakao, T. Hanashima, K. Munakata, T. Moyoshi, T. Kuroda, H. Kimura, T. Sakakura, C.H. Lee, M. Takahashi, K.I. Ohshima, T. Kiyotani, Y. Noda, M. Arai, J. Appl. Crystallogr. 49 (2016) 120. [6] K. Kusaka, T. Hosoya, I. Tanaka, N. Niimura, K. Kurihara, T. Ohhara, J. Phys.: Conf. Ser. 251 (2010) 012031. [7] R. Bolotovsky, M.A. White, A. Darovsky, P. Coppens, J. Appl. Crystallogr. 28 (1995) 86. [8] J. Peters, J. Appl. Crystallogr. 36 (2003) 1475. [9] A.C.Larsen, R.B.V.Dreele, GSAS – General Structure Analysis, System, Los Alamos National Laboratory, Report LAUR86-748 [10] D.A.Keen, C.C.Wilson, Rutherford Appleton Laboratory, Report RAL-TR-96-083

Acknowledgements I would like to thank Dr. Silvia Capelli for encouraging me to write this paper which has been long overdue. My thanks go to Dr. Devinder Sivia for carefully checking my calculations and his early efforts into implementing such a function. References [1] D.A. Keen, M.J. Gutmann, C.C. Wilson, J. Appl. Crystallogr. 39 (2006) 714. [2] G. Jogl, X. Wang, S.A. Mason, A. Kovalevsky, M. Mustyakimov, Z. Fisher, C. Hoffman, C. Kratky, P. Langan, Acta Crystallogr. Sect. D. 67 (2011) 584.


Nuclear Instruments and Methods in Physics Research A (xxxx) xxxx–xxxx

M.J. Gutmann

York, 1996. [15] C.L. Bull, M.W. Johnson, H. Hamidov, K. Komatsu, M. Guthrie, M.J. Gutmann, J.S. Loveday, R.J. Nelmes, J. Appl. Crystallogr. 47 (2014) 974. [16] P. Krooss, T. Niendorf, P.M. Kadletz, C. Somsen, M.J. Gutmann, Y.I. Chumlyakov, W.W. Schmahl, G. Eggeler, H.J. Maier, Shape Mem. Superelasticity 1 (2015) 6. [17] M.J. Gutmann, G. Graziano, S. Mukhopadhyay, K. Refson, M. von Zimmermann, J. Appl. Crystallogr. 48 (2015) 1122.

[11] K. Tomoyori, K. Kusaka, T. Yamada, T. Hosoya, T. Ohhara, K. Kurihara, I. Tanaka, M. Katagiri, N. Niimura, Nucl. Instrum. Methods Phys. Res. A 723 (2013) 128. [12] A.J. Schultz, M.R.V. Jørgensen, X. Wang, R.L. Mikkelson, D.J. Mikkelson, V.E. Lynch, P.F. Peterson, M.L. Green, C.M. Hoffmann, J. Appl. Crystallogr. 47 (2014) 915. [13] M.J. Gutmann, SXD2001, ISIS Facility, Rutherford Appleton Laboratory, Oxfordshire, England, 2005. [14] D.S. Sivia, Data Analysis – A Bayesian Tutorial, Oxford University Press Inc., New