Rietveld refinement is a technique described by Hugo Rietveld for use in the characterisation of crystalline materials. The neutron and X-ray diffraction of powder samples results in a pattern characterised by reflections (peaks in intensity) at certain positions. The height, width and position of these reflections can be used to determine many aspects of the material's structure.
The Rietveld method uses a least squares approach to refine a theoretical line profile until it
matches the measured profile. The introduction of this technique was a significant step forward in the
diffraction analysis of powder samples as, unlike other techniques at that time, it was able to deal reliably with strongly overlapping reflections.
The method was first implemented in 1967, and reported in 1969 for the diffraction of monochromatic neutrons where the reflection-position is reported in terms of the Bragg angle, 2θ. This terminology will be used here although the technique is
equally applicable to alternative scales such as x-ray energy or neutron time-of-flight. The only wavelength and technique independent scale is in reciprocal space units or momentum transfer Q, which is historically rarely used in powder diffraction but very common in all other diffraction and optics techniques. The relation is
:<math> Q = \frac {4 \pi \sin \left ( \theta \right )}{\lambda}.</math>
Introduction
The most common powder X-ray diffraction (XRD) refinement technique used today is based on the method proposed in the 1960s by Hugo Rietveld.
Powder diffraction profiles: peak positions and shapes
Before exploring Rietveld refinement, it is necessary to establish a greater understanding of powder diffraction data and what information is encoded therein in order to establish a notion of how to create a model of a diffraction pattern, which is of course necessary in Rietveld refinement. A typical diffraction pattern can be described by the positions, shapes, and intensities of multiple Bragg reflections. Each of the three mentioned properties encodes some information relating to the crystal structure, the properties of the sample, and the properties of the instrumentation. Some of these contributions are shown in Table 1, below.
<br />
{| class="wikitable"
|+Powder diffraction pattern as a function of various crystal structure, specimen, and instrumental parameters
!Pattern component
!Crystal structure
!Specimen property
!Instrumental parameter
|-
|Peak position
|Unit cell parameters
(a, b, c, α, β, γ)
|
- Absorption
- Porosity
|
- Radiation (wavelength),
- Instrument/sample alignment
- Axial divergence of the beam
|-
|Peak intensity
|Atomic parameters
(x, y, z, B, etc.)
|
- Preferred orientation
- Absorption
- Porosity
|
- Geometry and configuration
- Radiation (Lorentz polarization)
|-
|Peak shape
|
- Crystallinity
- Disorder
- Defects
|
- Grain size
- Strain
- Stress
|
- Radiation (spectral purity)
- Geometry
- Beam Conditioning
|}
The structure of a powder pattern is essentially defined by instrumental parameters and two crystallographic parameters: unit cell dimensions, and atomic content and coordination. So, a powder pattern model can be constructed as follows:
- Establish peak positions: Bragg peak positions are established from Bragg's law using the wavelength and d-spacing for a given unit cell.
- Determine peak intensity: Intensity depends on the structure factor, and can be calculated from the structural model for individual peaks. This requires knowledge of the specific atomic coordination in the unit cell and geometrical parameters.
- Peak shape for individual Bragg peaks: Represented by functions of the FWHM (which vary with Bragg angle) called the peak shape functions. Realistically ab initio modelling is difficult, and so empirically selected peak shape functions and parameters are used for modelling.
- Sum: The individual peak shape functions are summed and added to a background function, leaving behind the resultant powder pattern.
It is easy to model a powder pattern given the crystal structure of a material. The opposite, determining the crystal structure from a powder pattern, is much more complicated. A brief explanation of the process follows, though it is not the focus of this article.
To determine structure from a powder diffraction pattern the following steps should be taken. First, Bragg peak positions and intensities should be found by fitting to a peak shape function including background. Next, peak positions should be indexed and used to determine unit cell parameters, symmetry, and content. Third, peak intensities determine space group symmetry and atomic coordination. Finally, the model is used to refine all crystallographic and peak shape function parameters. To do this successfully, there is a requirement for excellent data which means good resolution, low background, and a large angular range.
Peak shape functions
For general application of the Rietveld method, irrespective of the software used, the observed Bragg peaks in a powder diffraction pattern are best described by the so-called peak shape function (PSF). The PSF is a convolution of three functions: the instrumental broadening <math>\Omega(\theta)</math>, wavelength dispersion <math>\Lambda(\theta)</math>, and the specimen function <math>\Psi(\theta)</math>, with the addition of a background function, <math>b(\theta)</math>. It is represented as follows:
:<math>PSF(\theta) = \Omega(\theta)\otimes\Lambda(\theta)\otimes\Psi(\theta) + b(\theta)</math>,
where <math>\otimes</math> denotes convolution, which is defined for two functions <math>f</math> and <math>g</math> as an integral:
:<math>f(t)\otimes g(t) = \int_{-\infty}^{\infty} f(\tau)g(t-\tau)d\tau = \int_{-\infty}^{\infty} g(\tau)f(t-\tau)d\tau </math>
The instrumental function depends on the location and geometry of the source, monochromator, and sample. Wavelength function accounts for the distribution of the wavelengths in the source, and varies with the nature of the source and monochromatizing technique. The specimen function depends on several things. First is dynamic scattering, and secondly the physical properties of the sample such as crystallite size, and microstrain.
A short aside: unlike the other contributions, those from the specimen function can be interesting in materials characterization. As such, average crystallite size, <math>\tau</math>, and microstrain, <math>\varepsilon</math>, effects on Bragg peak broadening, <math>\beta</math> (in radians), can be described as follows, where <math>k</math> is a constant:
:<math>\beta = \frac{\lambda}{\tau \cdot \cos\theta} </math> and <math> \beta = \kappa \cdot \varepsilon \cdot \tan\theta </math>.
Returning to the peak shape function, the goal is to correctly model the Bragg peaks which exist in the observed powder diffraction data. In the most general form, the intensity, <math>Y(i)</math>, of the <math>i^\text{th} </math> point (<math>1\leq i \leq n</math>, where <math>n</math> is the number of measured points) is the sum of the contributions <math>y_k</math> from the <math>m</math> overlapped Bragg peaks (<math>1 \leq k \leq m</math>), and the background, <math>b(i)</math>, and can be described as follows:
:<math>Y(i) = b(i)+ \sum_{k=1}^{m}{I_k[y_k(x_k)]} </math>
where <math> I_k </math> is the intensity of the <math> k^\text{th} </math> Bragg peak, and <math> x_i = 2\theta_i - 2\theta_k </math>. Since <math> I_k </math> is a multiplier, it is possible to analyze the behaviour of different normalized peak functions <math> y(x) </math> independently of peak intensity, under the condition that the integral over infinity of the PSF is unity. There are various functions that can be chosen to do this with varying degrees of complexity. The most basic functions used in this way to represent Bragg reflections are the Gauss, and Lorentzian functions. Most commonly though, is the pseudo-Voigt function, a weighted sum of the former two (the full Voigt profile is a convolution of the two, but is computationally more demanding). The pseudo-Voigt profile is the most common and is the basis for most other PSF's. The pseudo-Voigt function can be represented as:
:<math> y(x) = V_p(x) = n * G(x) + (1-n) * L(x) </math>,
where
:<math> G(x) = \frac{C_G^\frac{1}{2{\sqrt{\pi}H}e^{-C_Gx^2} </math>
and
:<math>L(x) = \frac{C_L^\frac{1}{2{\sqrt{\pi}H'} \left(1 + C_L x^2\right)^{-1}</math>
are the Gaussian and Lorentzian contributions, respectively.
Thus,
:<math> V_p(x) = \eta \frac{C_G^{\frac{1}{2}{\sqrt{\pi}H}e^{-C_Gx^2} + (1-\eta) \frac{C_L^\frac{1}{2{\sqrt{\pi}H'}(1+C_Lx^2)^{-1}. </math>
where:
- <math>H</math> and <math>H'</math> are the full widths at half maximum (FWHM)
- <math>x = \frac{2\theta_i - 2\theta_k}{H_k}</math> is essentially the Bragg angle of the <math>i^\text{th}</math> point in the powder pattern with its origin in the position of the <math>k^\text{th}</math> peak divided by the peak's FWHM.
- <math>C_G = 4\ln2</math>, <math> C_L = 4 </math> and <math display="inline"> \frac{C_G^{\frac{1}{2}{\sqrt{\pi}H} </math> and <math display="inline">\frac{C_L^{\frac{1}{2}{\sqrt{\pi}H'} </math> are normalization factors such that <math display="inline"> \int_{-\infty}^{\infty}G(x)dx = 1 </math> and <math display="inline"> \int_{-\infty}^{\infty}L(x)dx = 1 </math> respectively.
- <math> H^2 = U\tan^2\theta + V\tan\theta + W </math>, known as the Caglioti formula, is the FWHM as a function of <math> \theta </math> for Gauss, and pseudo-Voigt profiles. <math>U</math>, <math>V</math>, and <math>W</math> are free parameters.
- <math display="inline"> H' = \frac{X}{\cos\theta} + Y\tan\theta </math> is the FWHM vs. <math> 2\theta </math> for the Lorentz function. <math>X</math> and <math>Y</math> are free variables
- <math> \eta = \eta_0 + \eta_1 2\theta + \eta_2 \theta^2 </math>, where <math> 0 \leq \eta \leq 1 </math> is the pseudo-Voigt mixing parameter, and <math> \eta_{0,1,2} </math> are free variables.
The pseudo-Voigt function, like the Gaussian and Lorentz functions, is a centrosymmetric function, and as such does not model asymmetry. This can be problematic for non-ideal powder XRD data, such as those collected at synchrotron radiation sources, which generally exhibit asymmetry due to the use of multiple focusing optics.
The Finger–Cox–Jephcoat function is similar to the pseudo-Voigt, but has better handling of asymmetry, which is treated in terms of axial divergence. The function is a convolution of pseudo-Voigt with the intersection of the diffraction cone and a finite receiving slit length using two geometrical parameters, <math>S/L</math>, and <math>H/L</math>, where <math>S</math> and <math>H</math> are the sample and the detector slit dimensions in the direction parallel to the goniometer axis, and <math>L</math> is the goniometer radius.
Peak shape as described in Rietveld's paper
The shape of a powder diffraction reflection is influenced by the characteristics of the beam, the experimental
arrangement, and the sample size and shape. In the case of monochromatic neutron sources the convolution
of the various effects has been found to result in a reflex almost exactly Gaussian in shape. If this
distribution is assumed then the contribution of a given reflection to the profile <math>y_i</math> at position <math>2 \theta_i</math>
is:
:<math>
y_i = I_k \exp \left [ \frac{-4 \ln \left (2 \right )}{H_k^2} \left (2\theta_i - 2\theta_k \right )^2 \right ]
</math>
where <math>H_k</math> is the full width at half peak height (full-width half-maximum), <math>2\theta_k</math> is the center of the reflex, and <math>I_k</math> is the calculated intensity of the reflex (determined from the structure factor,
the Lorentz factor, and multiplicity of the reflection).
At very low diffraction angles the reflections may acquire an asymmetry due to the vertical divergence of the beam.
Rietveld used a semi-empirical correction factor, <math>A_s</math> to account for this asymmetry:
:<math> A_s = 1 - \left [ \frac {P \left (2\theta_i - 2\theta_k \right )^2}{\tan \theta_k} \right ] </math>
where <math>P</math> is the asymmetry factor and <math>s</math> is <math>+1</math>, <math>0</math>, or <math>-1</math> depending on the difference <math>2\theta_i-2\theta_k</math> being positive, zero, or negative respectively.
At a given position more than one diffraction peak may contribute to the profile. The intensity is simply the sum of all reflections contributing at the point <math>2\theta_i</math>.
Integrated intensity
For a Bragg peak <math>(hkl)</math>, the observed integrated intensity, <math>I_{hkl}</math>, as determined from numerical integration is
:<math>I_{hkl} = \sum_{i=1}^j(Y_i^{obs}-b_i)</math>,
where <math>j</math> is the total number of data points in the range of the Bragg peak. The integrated intensity depends on multiple factors, and can be expressed as the following product:
:<math>I_{hkl} = K \times p_{hkl} \times L_\theta \times P_\theta \times A_\theta \times T_{hkl} \times E_{hkl} \times |F_{hkl}|^2</math>
where:
- <math>K</math>: scale factor
- <math>p_{hkl}</math>: multiplicity factor, which accounts for symmetrically equivalent points in the reciprocal lattice
- <math>L_\theta</math>: Lorentz multiplier, defined by diffraction geometry
- <math>P_\theta</math>: polarization factor
- <math>A_\theta</math>: absorption multiplier
- <math>T_{hkl}</math>: preferred orientation factor
- <math>E_{hkl}</math>: extinction factor (often neglected as it is usually insignificant in powders)
- <math>F_{hkl}</math>: structure factor as determined by the crystal structure of the material
Peak width as described in Rietveld's paper
The width of the diffraction peaks are found to broaden at higher Bragg angles. This angular dependency was originally represented by
:<math> H_k^2 = U \tan^2 \theta_k + V \tan \theta_k + W </math>
where <math>U</math>, <math>V</math>, and <math>W</math> are the half-width parameters and may be refined during the fit.
Preferred orientation
In powder samples there is a tendency for plate- or rod-like crystallites to align themselves along the axis of a cylindrical sample holder. In solid polycrystalline samples the production of the material may result in greater volume fraction of certain crystal orientations (commonly referred to as texture). In such cases the reflex intensities will vary from that predicted for a completely random distribution. Rietveld allowed for moderate cases of the former by introducing a correction factor:
:<math> I_\text{corr} = I_\text{obs} \exp \left (-G\alpha^2 \right ) </math>
where <math>I_\text{obs}</math> is the intensity expected for a random sample, <math>G</math> is the preferred orientation parameter and <math>\alpha</math> is the acute angle between the scattering vector and the normal of the crystallites.
Refinement
The principle of the Rietveld method is to minimize a function <math>M</math> which analyzes the difference between a calculated profile <math>y^\text{calc}</math> and the observed data <math>y^\text{obs}</math>. Rietveld defined such an equation as:
:<math> M = \sum_{i} W_i \left \{ y_i^\text{obs} - \frac{1}{c} y_i^\text{calc} \right \}^2 </math>
where <math>W_i</math> is the statistical weight and <math>c</math> is an overall scale factor such that <math>y^\text{calc} = c y^\text{obs}</math>.
Least squares method
The fitting method used in Rietveld refinement is the non-linear least squares approach. A detailed derivation of non-linear least squares fitting will not be given here. Further detail can be found in Chapter 6 of Pecharsky and Zavalij's text 12 . There are a few things to note however. First, non-linear least squares fitting has an iterative nature for which convergence may be difficult to achieve if the initial approximation is too far from correct, or when the minimized function is poorly defined. The latter occurs when correlated parameters are being refined at the same time, which may result in divergence and instability of the minimization. This iterative nature also means that convergence to a solution does not occur immediately for the method is not exact. Each iteration depends on the results of the last which dictate the new set of parameters used for refinement. Thus, multiple refinement iterations are required to eventually converge to a possible solution.
Rietveld method basics
Using non-linear least squares minimization, the following system is solved:
:<math>\begin{pmatrix} Y_i^\text{calc} = kY_i^\text{obs} \\\vdots\\ Y_n^\text{calc} = kY_n^\text{obs} \end{pmatrix}</math>
where <math>Y_i^\text{calc}</math> is the calculated intensity and <math>Y_i^\text{obs}</math> is the observed intensity of a point <math>i</math> in the powder pattern, <math>k</math>, is a scale factor, and <math>n</math> is the number of measured data points. The minimized function is given by:
:<math>\Phi= \sum_{i=1}^nw_i(Y_i^\text{obs}-Y_i^\text{calc})^2</math>
where <math>w_i</math> is the weight, and <math>k</math> from the previous equation is unity (since <math>k</math> is usually absorbed in the phase scale factor). The summation extends to all <math>n</math> data points. Considering the peak shape functions and accounting for the overlapping of Bragg peaks because of the one-dimensionality of XRD data, the expanded form of the above equation for the case of a single phase measured with a single wavelength becomes:
:<math>\Phi= \sum_{i=1}^nw_i\biggl( Y_i^\text{obs}-\Bigl(b_i+K\sum_{j=1}^mI_jy_j(x_j)\Bigr)\biggr)^2</math>
where:
- <math>b_i</math> is the background at the <math>i^\text{th}</math> data point.
- <math>K</math> is the phase scale factor.
- <math>m</math> is the number of Bragg reflections contributing to the intensity of the <math>i^\text{th}</math> reflection.
- <math>I_j</math> is the integrated intensity of the <math>j^\text{th}</math> Bragg peak.
- <math>y_i(x_i)</math> is the peak shape function.
For a material that contains several phases (<math>p</math>), the contribution from each is accounted for by modifying the above equation as follows:
:<math>\Phi= \sum_{i=1}^nw_i\biggl( Y_i^\text{obs}-\Bigl(b_i+ \sum_{l=1}^pK_l\sum_{j=1}^mI_{l,j}y_{l,j}(x_{l,j})\Bigr)\biggr)^2</math>
It can easily be seen from the above equations that experimentally minimizing the background, which holds no useful structural information, is paramount for a successful profile fitting. For a low background, the functions are defined by contributions from the integrated intensities and peak shape parameters. But with a high background, the function being minimized depends on the adequacy of the background and not integrated intensities or peak shapes. Thus, a structure refinement cannot adequately yield structural information in the presence of a large background.
It is also worth noting the increased complexity brought forth by the presence of multiple phases. Each additional phase adds to the fitting, more Bragg peaks, and another scale factor tied to corresponding structural parameters, and peak shape. Mathematically they are easily accounted for, but practically, due to the finite accuracy and limited resolution of experimental data, each new phase can lower the quality and stability of the refinement. It is advantageous to use single phase materials when interested in finding precise structural parameters of a material. However, since the scale factors of each phase are determined independently, Rietveld refinement of multi phase materials can quantitatively examine the mixing ratio of each phase in the material.
Refinement parameters
Background
Generally, the background is calculated as a Chebyshev polynomial. In GSAS and GSAS-II they appear as follows. Again, background is treated as a Chebyshev polynomial of the first kind ("Handbook of Mathematical Functions", M. Abramowitz and IA. Stegun, Ch. 22), with intensity given by:
:<math>I_i = \sum_{j=1}^NP_jT'_{j-1}</math>
where <math>T'_{j-1}</math> are the coefficients of the Chebyshev polynomial taken from Table 22.3, pg. 795 of the Handbook. The coefficients have the form:
:<math>T'_n = \sum_{m=0}^{i-1}C_mX^m</math>
and the values for <math>C_m</math> are found in the Handbook. The angular range (<math>2\theta</math>) is converted to <math>X</math> to make the Chebyshev polynomial orthogonal by
:<math>X = \frac{2}{\tau}-1</math>
And, the orthogonal range for this function is –1 to +1.
Sample displacement: sample transparency and zero-shift corrections
Different correction factors have been developed to account for the specimen-detector displacement in Debye-Scherrer (transmission) and Bragg-Brentano (reflection) geometries. Correction factors have been published for flat tilted and non-tilted 2D detectors as well.
Traditional X-ray diffraction experiments place the detector normal to the incident X-ray beam, centering the detector on the beam or placing the incident beam at the bottom edge of the detector. In these experiments, a unique configuration is taken by aligning the detector such that the incident X-rays are still hitting the bottom-most edge of the detector. However, tilting the detector forward such that the topmost edge of it is directly above the sample provides a number of benefits such as higher quality and better resolved radial distribution functions than in typical traditional geometries, improved dynamic range, increased reciprocal space resolution per pixel for low angle scattering and access to higher reciprocal space values at lower energies. Tilting the detector does not have much impact when carrying out SAXS experiments, however it has a huge impact on WAXS measurements in large-scale facilities (synchrotrons).
thumb|Sample-detector displacement correction with flat 2D detectors and non-conventional geometries
The appropriate function for the correction of <math>2 \theta_c</math> peak positions due to sample displacement depends upon the geometry of the instrument. Guzmán’s equation The most popular and conventional figure of merit used is the goodness of fit which should approach unity given a perfect fit, though this is rarely the case. In practice, the best way to assess quality is a visual analysis of the fit by plotting the difference between the observed and calculated data plotted on the same scale.