An Implementation of the Fundamental Parameters Approach for Analysis of X-ray Powder Diffraction Line Profiles

This work presents an open implementation of the Fundamental Parameters Approach (FPA) models for analysis of X-ray powder diffraction line profiles. The original literature describing these models was examined and code was developed to allow for their use within a Python based least squares refinement algorithm. The NIST interest in the FPA method is specific to its ability to account for the optical aberrations of the powder diffraction experiment allowing for an accurate assessment of lattice parameter values. Lattice parameters are one of the primary certified measurands of NIST Standard Reference Materials (SRMs) for powder diffraction. Lattice parameter values obtained from analysis of data from SRMs 640e and 660c using both the NIST FPA Python code and the proprietary, commercial code Topas, that constitutes the only other actively supported, complete implementation of FPA models within a least-squares data analysis environment, agreed to within 2 fm. This level of agreement demonstrates that both the NIST code and Topas constitute an accurate implementation of published FPA models.


Introduction
Measured X-ray powder diffraction line profiles are affected by the geometry of the diffractometer, the shape of the X-ray emission spectrum and the physical characteristics of the sample as shown diagrammatically in Fig. 1. Certain aberrations embodied in the geometric instrument profile are highly asymmetric and will displace the observed position of profile maxima; therefore, these maxima are not indicative of the true lattice spacing. Furthermore, the level and direction of the displacement can vary dramatically as a function of diffraction angle. High-accuracy determination of the crystallite lattice spacing from measured line profiles requires accounting for such effects in the model for line profile shape. The Fundamental Parameters Approach (FPA) for the description of X-ray line profiles is a convolutionbased line profile modeling method that describes the measured line profiles as the convolution of peak profiles representing the emission spectrum with a number of aberration functions, each representing a certain aspect of the instrument configuration or sample microstructure that affects the measured peak position and/or shape. The parameters of FPA models are entirely interpretable in terms of the physical  [ 1] ε −  is the last element [ : ] i j ε  the elements of an array ε  with indices between i and j−1 (inclusive), so it has length j−i x ε +  an operation between an array and a scalar operates element-by-element on the array with the scalar x ε +   an operation between two arrays is done element-by-element ( ) f ε  a function applied to an array is an array of the same length with the function applied to each element # in pseudo-code sections, everything after this on a line is a comment scaling we will present all equations below in a manner that is mostly compatible with the usage established by Topas. Lorentzian widths and Gaussian widths are expressed as the full-width at half-maximum (FWHM) of the peak shape. However, all lengths are uniformly scaled; any consistent unit of length can be used, but all lengths must be the same units. The reference code we provide takes all angles in degrees, and converts them internally to radians.

Initialization of Parameters
To start a calculation, we assume that the result will be a peak shape, uniformly gridded in 2θ space, centered at 2θ 0 , with N 2θ points (which will typically be even and either a power of two or, for modern Fast Fourier Transforms (FFT) packages, products of powers of 2, 3, 5, and 7, typically), and with full width W. To match the standard of typical Fourier transform package implementations, using a common convention for transforms of purely real data sets, we will carry only the coefficients corresponding to non-negative angular frequencies. We will assume the length of the ω  array is 2 /2 1 N θ + , in which the highest frequency cosine component is put at the end, rather than folded back into the zero-frequency bin as is the result of full, complex FFTs. The ω  array is initialized, then, to 226 http://dx.doi.org/10.6028/jres.120.014 where W is in radians.

Source Spectrum and Crystal Size
We handle the convolutions due to the emission spectrum and due to the crystal size parameters first, and together, because it is numerically efficient to do so. We initialize the Fourier transform buffer F  that will be multiplied with all the other convolutions with the sum of the transforms of the emission spectrum broadened by the peak size. These two are kept together because they each have Lorentzian and Gaussian components that are easily combined. It is important to note that these crystallite size contributions have nothing to do with the actual emission spectrum; the grouping is purely historical and convenient.
For an emission spectrum which is a set of n λ lines, indexed by k, each with Lorentzian FWHM l k , Gaussian FWHM g k , intensity I k , and wavelength k λ , and a Gaussian crystal size component S G and a Lorentzian crystal size component S L , and a reference wavelength 0 λ which is typically that of the strongest line in the spectrum, we define: Then the transform buffer F  is initialized via: Care must be exercised with respect to the sign on the complex term in the exponential. Some Fourier transform packages define this differently. It is suggested that if a code is written, the order of the lines in the spectrum be verified. If the spectrum is backwards in the diffraction pattern, this sign should be switched.

Axial Divergence
The axial divergence model of Cheary and Coelho [5] (referred to as the "full" axial divergence model in Topas) is effective for the computation of this complicated function for a wide variety of situations. However, it is very intricate, involving many tests for various boundary conditions, and if implemented without attention to numerical issues can be computationally expensive. We present an implementation that is functionally equivalent, but has many of the boundary conditions clarified, and which is less subject to numerical issues than a literal implementation would be. In particular, we handle the 1/ x divergences so that the function can be evaluated more coarsely without loss of smoothness.
For the purposes of a practical implementation of this model, we need to be able to calculate 2 ( , ) I β ε per Sec. 4.2 of [5] and then, ultimately, the full intensity profile including the effects of Soller slits via integration of their equation 27. The remainder of this section will present this in detail, concentrating on a clean algorithmic implementation, with all theoretical details being referred back to [5,8]. The first step of the calculation is to set up some parameters that will be needed (with locations in [5] where appropriate): Then, following eqs. 15a, b, c, d of [5], we compute the parameters 0 Z + and 0 Z − conditionally on β : Then, from eqs. 18a, b of [5] (with corrections to 18b),  . The problem is then divided into two major domains, each with three minor ones. This is most easily done following Table 1 of [5]. This has been amended to include omitted reflections for formulas when 2 > 90 θ  . We present this in two formats: a pseudocode block in algorithm 1, and Table 1.

Computing Components of Axial Divergence Shape for Fixed β
Using the parameters from algorithm 1 or Table 1, we need to set up the equations from Table 1 of [5] for F 1 , F 2 , F 3 , and F 4 . These will be used exactly as defined in the original work, with r β defining the β range to select the equations. However, this is the component of the computation where the most critical numerical issues must be addressed. First, we assume that the angle offset ε is defined on a uniform grid centered at 0 and running from w − to w where w is the half-width of a window on which the axial divergence function in being computed. It will be treated as an array ε  and stored as an array of n elements. All of the F functions compute pieces of a function 0 0 / | | y k ε ε + − . This computation may well include the endpoint where 0 = ε ε ,where this function diverges. Also, due to the discrete binning of ε , it may include a point at which the argument (without the absolute value) is in fact negative, but should be truncated to 0. Further, due to the discrete binning, the sum of all of the sampled values of the function may not add up to the integral of the function over the bounds, resulting in inaccuracies in the total X-ray peak intensity, especially in the http://dx.doi.org/10.6028/jres.120.014 Table 1. Selection of computation boundaries and β ranges

Condition 1 and
case of sampling on a relatively coarse grid for computational speed. Finally, again due to discrete sampling, the first moment (centroid) of the distribution may not be exactly correct. This is particularly critical, since shifts in this result in inaccuracy of peak positions. Such shifts can be reduced by using finer computational grids, but this approach is very inefficient. Therefore, we adjust the result so that, in all cases, it has the exact centroid one would expect in the continuum limit. All of these issues can be addressed with a single 'helper' function F 0 which will be used to construct F 1,2,3,4 in a unified manner. This is presented in the next section.

Helper Function F 0
This function will take as formal arguments describe it algorithmically, the details of it are in the Python implementation of the FPA provided as ancillary materials with this paper. The function 'axial_helper' in the python code is the implementation of F 0 .

Computing Complete Axial Divergence for Fixed β
With the assistance of the function in 2.4.2, the bookkeeping in Table 1 of [5] is straightforward to implement. The procedure for doing this is shown as pseudo-code in 2.4.4. When this process is complete, one has in hand the arrays representing the functions 2 I + and 2 I − for a given β and 0 2θ . Table 1 Computation for I 2 We define the various functions needed for Table 1 in terms of the helper F 0 . Note that F 0 modifies the dst   argument in place, and returns the lower and upper bounds of the part of the array which is modified.

Carrying Out
or, equivalently: Then, we create arrays for 2 I + and 2 I − which will be the accumulators as defined in caption of Table 1 of [5]. Also, we create an empty list which will accumulate index bounds. The semantics of this list will be that the + = operator concatenates the elements returned by the function onto the end. The algorithm is shown in algorithm 2.

Computing I 3 by Integrating I 2
The previous two sections have presented the most complex part of the computation of the axial divergence peak shape. This section presents an implementation of Sec. 5 of [5], the computation of 3 ( , ) I ε β and the integral of it over all allowed β to get the total line shape. First, we create functions representing the transmission of Soller slits of angular full-width β max for the incident Soller slit and angular full-width max γ for the receiving Soller slit. These are from equation 24a, b of [5], slightly rewritten: Algorithm 3 shows how to carry out the integral over β.

Applying the Axial Divergence Convolution
To apply this convolution, which has been generated in real space, it must be transformed to Fourier space: Because of the way the peak is centered, this transform has an alternating sign of 1 − across its values.
All odd-numbered bins of ( ) g ω  will be multiplied by −1. Then, F  will be multiplied by ( ) g ω  .

Sample Offset and Zero Angle
The correction for the sample surface not lying in the equatorial plane of the diffractometer can be combined with the shift due to the zero angle error. If the sample is offset by 0 z , and the zero angle offset and the convolver is: and F  is multiplied by ( ) g ω  . http://dx.doi.org/10.6028/jres.120.014

Source Width and Tube Tails
The correction for the broad shoulders on the spatial distribution of emission for fine-focus X-ray tubes can be computed according to [7] and Sec. 4.1 of [8]. Define m w as the width of the central peak, l w as the distance to the low-angle side of the shoulder, h w as the distance to the high-angle side, and t I as the intensity of the tails, then, and the convolution is and F  will be multiplied by this ( ) g ω  .

Receiver Slit Equatorial Height
A receiver slit of full height r h creates a rectangular convolution of angular full width / r h R (where R is the diffractometer radius). From 42, the Fourier space representation is: which will be multiplied into F  .

Flat Specimen/Equatorial Divergence Slit Size
From eqs. 9 and 10 of [8], the correction for the flat specimen error is where α is the equatorial divergence angle of the X-ray beam. Then, the convolution function is This is of the same form as 0 F , used in the axial divergence calculation (see Sec. 2.4.2), so we will use that function. Note that this is a computation in real space. This is shown in algorithm 4.
Because of the way the peak is centered, this transform has an alternating sign of 1 − across its values. All odd-numbered bins of ( ) g ω  will be multiplied by −1. Then, F  will be multiplied by ( ) g ω  .

Specimen Transparency
From equation 12 of [8], with a correction to δ , the convolution due to the finite interaction depth in the target is: where m is the absorption coefficient of the sample, measured in units consistent with R, so if R is in mm, m would be in mm −1 , and T is the sample thickness in the same units as R. From 48, the Fourier transform of this expression is: Then, F  will be multiplied by ( ) g ω  .

θ Ω ≠ )
If the specimen angle Ω is offset from θ , a defocusing correction appears, per equation 15 of [8]. It is a rectangular convolution of angular full width ( ) where α is the equatorial divergence angle. For small θ − Ω , this also reduces to: From 42, the Fourier space representation is: which will be multiplied into F  .

Silicon Position Sensitive Detector (Si PSD)
A Si PSD, which has a finite window width, suffers defocusing due to two effects. The first is the effect discussed above, due to the sample angle Ω being different from the diffraction angle θ , resulting in a violation of the expected parafocusing optics. The second effect is due to the flat face of the detector itself; active pixels are not located on the radius defined by the diffractometer configuration. Unlike corrections for older PSDs [15], which include parallax due to the long absorption length in a gas-filled detector, we assume the detector is planar and has no effective depth.
The exact expression for a ray starting at the source at an angle α from the center ray, being diffracted by an angle 2θ , and striking the detector face which is centered at 2θ ε + , is a ray which intersects the detector face at a position y : This can be expanded as a series in ε and α , which is: The first term is just the expected offset of the peak on the detector face, and is not an aberration. The second term contains two components, the first of which is exactly that of 23. The second component depends on 2 ε . In the limit that θ → Ω , it reduces to simply 2 /2 ε . Thus, to the width for a high-angle peak. We present below an exact solution to the case in which the 2 ε is ignored. In the case it is needed, the integral below needs to be carried out numerically. If the detector window being analyzed extends from distance 1 y to distance 2 y above the centerline, and is centered so that 0 = θ Ω , where psd δ comes from 29. Note that the defocusing correction of Sec. 2.5.6 is symmetrical in 2θ , so it is only necessary to integrate over one side of the detector. The integral can be carried out analytically using the approximation of 24, resulting in an expression involving the sine integral function Si( ) Note that evaluation of this function involves a 0/0 where = 0 ω . However, the limit can easily be taken analytically. The zero bin of ( ) g ω  should be set to 2 1 ( )/ y y R − , which is just the angular length of the exposed detector face. Also note that the subtracted term in 32 vanishes if 1 = 0 y .
For a symmetrical detector, 1 = 0 y . However, some data analyses may split the data from the detector into a pattern using the central pixels of the detector to get a high-resolution result, and then use the remaining pixels to create another pattern which has lower angular resolution, but takes advantage of the counting statistics available over the full detector aperture. In this case, the 'central' pattern would use 1 = 0 y and 2 = hr y W where hr W is the half-height of the region to be sampled for high-resolution data. The 'outer' pattern would use 1 = hr y W and 2 = det y W where det W is the half-height of the full detector aperture. Such a split would permit nearly optimal use of the characteristics of the PSD, with high resolution and high active area.
The form of 32 gives direct guidance as to where the cutoff hr W should be made. For small arguments, Si( ) the Si PSD will not significantly broaden the overall response of the system as a result of this window selection. A more advanced approach would be to adjust this window width as a function of position in 2θ , so that as the peaks are dispersion-broadened at high angle, one would automatically use a wider window to collect more counts with no loss of resolution.
As with all the other corrections, F  will be multiplied by the resulting ( ) Aside on real-space solution If one is working in real space, rather than Fourier space, the integral of a top-hat function over widths DR δ of 24 can be carried out analytically to get the convolution function. As before, assuming the window extends from 1 y to 2, y , the result is: This expression includes a weak singularity at = 0 ε This can then have the forward difference computed, as in Sec. 2.4.2, to get the appropriate singularity-free binned function. A full code example is not provided, since it is essentially identical to the helper code.

Conversion of Transform Results to 2θ Space
In the previous sections, we have enumerated many convolutions that get accumulated multiplicatively into F  . Any other aberrations which are to be included can be done so in a similar manner. When everything is included, one needs to un-transform F  into real space to get , the nearly-final aberration function. We have been assuming the user has a Fourier transform library which provides a pair of properly-matched transform functions.
( )  takes a real array x  of length N and transforms it into the positive ω components of the Fourier transform, a complex array of length /2 1 N + . We apply its inverse, to get which takes /2 1 N + complex numbers and converts them back to N real elements. The one drawback to working in Fourier space is that all functions have built-in an implied periodicity of the length over which they are sampled. For many functions, this is not an issue, since they go quickly to zero near their boundaries, and so no wrap-around occurs. Unfortunately for the situation here, the http://dx.doi.org/10.6028/jres.120.014 aberration functions have extremely long Lorentzian tails, due to the contributions from the h k from 1. These tails always wrap around, which results in their inaccurate computation unless one spends much extra computational effort by computing the function over a very large interval in 2θ space. This problem has been addressed previously in a paper on the computation of Voigt functions by Fourier transform methods [16]. In short, the periodicity is equivalent to having computed the transform of an infinitely repeating comb of the line shape. Since the long tail is almost perfectly Lorentzian, one can subtract the infinite sum of Lorentzians from the computed aberration, which corrects the tail. This can be done in closed form. From equation 7 of [16], the sum is: where ∆ is the full width of the 2θ window, m is the computed centroid of the peak, and α is the halfwidth of the widest component of the Lorentzian. This function is normalized to unit area. Then, the corrected shape is: where A is the area of ( ) This makes a very good correction of the tails, assuming that the peak is not so asymmetrical that it has quite different amplitudes at the boundaries of the 2θ window. An example of the correction is shown in Fig. 2.   Fig. 2. Correction due to periodic Fourier transform, shown at low angle where the peak is very asymmetrical, and at mid-angle where it is nearly symmetrical. Note that for the left-hand case, the 2θ window is barely wide enough, so the peak tails are still very asymmetrical. http://dx.doi.org/10.6028/jres.120.014

Numerical Comparisons
The data in Secs. 3.1 and 3.2 are provided to allow one to compare implementations of the NIST FPAPC to the results from Topas. Section 3.3, with comparisons to data, is provided as general validation of the FPA for some important test cases. While the FPAPC is not a Rietveld code, in that it does not utilize a full structural mode, it can utilize space group symmetry (of various SRM materials) to constrain peak positions to a single lattice parameter, equivalent to a Pawley fit.

Simplified Source Spectrum, Variable Soller Slits
This section presents numerical and graphical comparisons of the output of FPAPC with that of Topas. The setup we are comparing is that of a diffractometer with the parameters shown in Table 2, using the point detector with the "full" model for axial divergence. These examples have the source spectrum artificially restricted, to make the effect of axial divergence and Soller slits more evident. We compute synthetic peak patterns for material with the characteristics of SRM 660c LaB 6 . Each data set is computed with different Soller slit settings, ranging from 2.5° full-width to 10.6° full-width. Tables 3,  4

Realistic Source Spectrum and Parameters
We now compute a sample with a more realistic, full spectrum as determined from fits to data from the NIST Johansson incident beam monochromator, as described in [1]. All parameters shown in Table 6 were fit by Topas and FPAPC. Table 7 shows the detailed errors, and Fig. 6 shows the peak shapes.

Comparison to Measurements
The most critical metric for comparison of the two programs is that of refined lattice parameter. SRMs 640e [17] and 660c [18] were certified in March, 2015 using data from the instrument described in [1]. The certification procedure involved the collection of twenty high-quality (24 hour scans) data sets for each of the two SRMs. These were analyzed independently using Topas with a Pawley analysis. With the NIST FPAPC, the twenty data sets were analyzed with a single, global refinement: specimen specific parameters, such as specimen displacement, were refined independently while instrument specific parameters, common across all data sets, were refined as single parameters. In Fig. 7 we show a typical fit to the data. Close correspondence between the instrument specific parameters obtained from Topas and those from FPAPC was observed. Additional testing indicated the residual error terms were not increased significantly by the variation of the instrument specific terms within the "window" of refined values obtained with the two codes. Instrument parameters, common to all data sets were then fixed at values that largely constituted the average values obtained with the two codes. This being done, the average of the lattice parameters values obtained from the average of the 20 independent Topas analyses, for both SRMs 640e and 660c, agreed with the corresponding global FPAPC values to within 2 fm.
Testing of the FPA model itself can be performed with an analysis of the variation in lattice parameter with reflection position in 2θ as reported previously [19]. This is a very sensitive test of the success of the FPA model as all the reflections in a pattern should give the same lattice parameter. Lattice parameter is the only property that is absolutely conserved across the entire pattern while profile asymmetry can vary dramatically in both degree and direction with 2θ. Again we used SRMs 640e and 660c for this purpose. In Fig. 7 we show the comparison of the peak positions from a globally refined lattice parameter, determined with FPAPC, with peak positions when refined independently. For a wide range of angles, from roughly 40° to 140° in 2θ, the corrections provided by the FPA are very good. There is, however, a clear, systematic tendency at low and high angles, where the peaks are most asymmetric, for the result to be biased. This is http://dx.doi.org/10.6028/jres.120.014 not understood, and is a matter of intense focus by the authors at this time. It is worth noting that the information of highest quality about the lattice parameter of the material comes from the peaks in the 50°-120° 2θ range, where contribution of aberrations are minimal and the angle is high enough that the contribution of a small angular error to the lattice parameter is minimized.
(a) Example fit of SRM660c data using FPAPC, with residuals (offset vertically for clarity).  SRM 1979 is being certified for the measurement of crystallite size from an analysis of profile broadening. It was prepared by decomposing zinc oxalate in a large-scale, NIST-built vacuum furnace using a heating schedule derived from the procedures outlined in [20,21]. The ZnO was then annealed in air to obtain two powders, one with an approximate crystallite size of 15 nm and a second one of 60 nm. In Fig. 8, we show a result of applying FPAPC code, extended to carry out the Scardi and Leoni model for log-normal crystallite size distributions [22] and the stacking fault density model of Warren [23]. This demonstrates that the algorithms described above can be extended to include complex models for material microstructure, many of which have natural representations in Fourier space. The breadth of the peaks in diffraction patterns from these ZnO materials varies widely due to the both crystallite size and the hkl specific stacking-faults.

Discussion
The Fundamental Parameters Approach to X-ray powder diffraction line profile analysis has played a central role in NIST powder diffraction SRM development since its inception with the aforementioned work of Cheary, Coelho and Cline (and collaborators on various SRM projects). We have demonstrated that the refined lattice parameter values obtained with our independently written NIST FPAPC and those of the commercial code Topas, that NIST has used since the year 2000 for SRM certification, agree to within 2 fm. This observation would confirm that both the NIST FPAPC and Topas are preforming in accordance to published FPA models. This conclusion is further supported by the data presented in Secs. 3.1 and 3.2 that illustrate that the form of the FPA profiles from the two programs are essentially identical. The equivalence of results between the NIST open implementation of FPA models and those from Topas enhances the transparency of the analyses performed in the certification of NIST SRMs for powder diffraction.
then its Fourier transform is:

One-Sided Exponential Tail
If ( ) = exp( ) f x x δ δ for 0 < < 0 Note that we present this in not-quite-normalized form. If 0 x → −∞ , it has unit area. We will use this in the case of finite sample thickness, in which case it should not be normalized to unit area, but just as shown here.

Translation of a Function
The Fourier transform of ( ) f x a + for any f is: Python name name in paper and explanation axial_helper F 0 ( ) _fxbuf no name in paper, a scratch buffer saved for memory efficiency searchsorted( ) a function which takes an ordered array and a list of values, and returns the bin index from the array for the number exactly matching each of the values, or returns the index of the bin to the right if the value falls between to bins.

Python Array Indexing
The most complex part of an efficient implementation of the FPA is the assignment of data to specific bins in a discretized representation of a computed diffraction profile. This is the reason for the complexity of F 0 . Because every computer language has its own conventions for how to index an array, and this is central to this work, we include here short notes about the Python language array indexing, to assist in translation to other languages. Thus, Python arrays: • are indexed in a zero-based manner; x[0] is the first element of an array • permit negative indexing to index from the last element; x[-1] is the last element of an array • can have a sub-array taken from them by indexing x [1:5]; this means the four elements of the array x[1], x[2], x [3], and x [4]. The length of a slice is the difference between the ending and starting indices. This is a feature of python indexing which can cause confusion, but is quite handy. This rule also implies that x[1:-1] indexes everything except the first and last elements of the array. • can have every second (e.g.) element, starting from 1 and ending before 9, indexed as x[1:9:2]; this extracts x[1], x [3], x [5], and x [7]. Note that x [9] is excluded since the upper end is never included. Thus, to negate all the odd-numbered elements of an array, x[1::2]*=-1. Note that empty index slots mean all possible, so this includes (potentially) the last element of the array. This is used extensively to re-phase Fourier transforms so that the zero position of the inverse transform is in the center of an array, rather than at the left edge, with negative positions wrapped to the right edge.

General Notes
The code, as provided in the supplementary material, is quite a bit more complex than a stripped-down reference implementation normally would be. This is due to the fact this version includes quite a bit of capability related to caching results which are likely to be re-used in the process of least-squares fitting of peaks. It is advised that the reader wishing to understand the underlying theory pay primary attention to the contents of the "conv_xxx" functions, which contain the machinery used to generate the convolutions, and to the various functions beginning with "axial_" which handle the real-space parts of the axial divergence integral.