Predicting XAFS scattering path cumulants and XAFS spectra for metals (Cu, Ni, Fe, Ti, Au) using molecular dynamics simulations

The ability of molecular dynamics (MD) simulations to support the analysis of X-ray absorption fine-structure (XAFS) data for metals is evaluated. The low-order cumulants (ΔR, σ², C₃) for XAFS scattering paths are calculated for the metals Cu, Ni, Fe, Ti and Au at 300 K using 28 interatomic potentials of the embedded-atom method type. The MD cumulant predictions were evaluated within a cumulant expansion XAFS fitting model, using global (path-independent) scaling factors. Direct simulations of the corresponding XAFS spectra, χ(R), are also performed using MD configurational data in combination with the FEFF ab initio code. The cumulant scaling parameters compensate for differences between the real and effective scattering path distributions, and for any errors that might exist in the MD predictions and in the experimental data. The fitted value of ΔR is susceptible to experimental errors and inadvertent lattice thermal expansion in the simulation crystallites. The unadjusted predictions of σ² vary in accuracy, but do not show a consistent bias for any metal except Au, for which all potentials overestimate σ². The unadjusted C₃ predictions produced by different potentials display only order-of-magnitude consistency. The accuracy of direct simulations of χ(R) for a given metal varies among the different potentials. For each of the metals Cu, Ni, Fe and Ti, one or more of the tested potentials was found to provide a reasonable simulation of χ(R). However, none of the potentials tested for Au was sufficiently accurate for this purpose.

1. Introduction

X-ray absorption fine-structure (XAFS) spectroscopy is a versatile technique for determination of local atomic structure (Bunker, 2010). Atomistic [Monte Carlo and molecular dynamics (MD)] simulations are increasingly used to calculate the input data (cumulants, synthetic XAFS spectra) required for XAFS structural analysis of solids (Di Cicco et al., 2002; a Beccara et al., 2003; Okamoto, 2004; Witkowska et al., 2006; Kuzmin & Evarestov, 2009; Kalinko et al., 2009; Higginbotham et al., 2009; Roscioni et al., 2011; Price et al., 2012). The primary motivation for using MD in XAFS analysis is often to generate trial structures that can be tested or refined during the course of a fitting procedure. Conversely, the theoretical basis of XAFS is sufficiently well established that XAFS can be used to test the structural and dynamical predictions of atomistic simulations (Hayes & Boyce, 1980; Mousseau & Thorpe, 1992; Newville, 1995; Edwards et al., 1997; Binsted et al., 2005). Further integration between XAFS and atomistic simulation techniques is desirable for both purposes. In particular, emerging materials technologies (Marletta et al., 2011; Hellborg et al., 2010) have led to an increasing need to characterize, through both experiment and simulation, the structures of nanoscale materials that are fabricated by modification of the near-surface region of solids (e.g. by ion bombardment or overlayer deposition) (MoberlyChan et al., 2007; Krasheninnikov & Nordlund, 2010; Baglin & Ila, 2011). This study benchmarks the accuracy of MD simulations for metallic solids against XAFS data, using a simulation methodology that would be appropriate for nanoscale materials applications.

The physics of a MD simulation is largely bound up in the properties of the interatomic potential. MD simulations of surface-terminated crystallites of the metals Cu, Ni, Fe, Ti and Au at 300 K are reported here for a representative selection of many-body interatomic potentials of the embedded atom method (EAM) type (Daw et al., 1993) that differ in terms of their functional forms and ranges. These metals include examples of the face-centered cubic (f.c.c.), body-centred cubic (b.c.c.) and hexagonal close packing (h.c.p.) crystal types. The structural properties used for fitting the EAM potentials are predominantly those of ideal crystals, whereas most nanoscale applications of XAFS (clusters, films, radiation effects) involve materials with defects or reduced dimensionality. The MD simulations were performed using surface-terminated crystallites, in order to better resemble typical experimental substrates and conditions. The presence of a surface is a necessary feature in any simulation that seeks to model materials modification by particle bombardment or deposition. The chosen boundary conditions provide a robust test of the ability of the EAM potentials to reproduce bulk structural properties in simulation crystallites that lack the rigid volume and shape constraints imposed by three-dimensional periodicity.

The MD configurational data that arise from the various metal/potential combinations are employed for calculation of (i) the cumulants of XAFS scattering path distributions, and (ii) synthetic XAFS spectra. The ability of MD to support the analysis of experimental XAFS data for these metals is evaluated for each potential by comparing the MD predictions of scattering path cumulants and synthetic XAFS spectra, respectively, to the corresponding data derived from experimental measurements. The potentials for a given metal are found to vary significantly in their ability to model XAFS data at the level of precision required for structural analysis.

2. Experimental

Transmission XAFS data for Cu, Ni, Fe (K-edge) and Au (L₃-edge) foil standards were recorded at the Pacific Northwest Consortium Collaborative Access Team (PNC-CAT) bending-magnet beamline (sector 20-BM) at the Advanced Photon Source (APS). Transmission XAFS K-edge data were also collected for Cu and Ti foil standards at the Australian National Beamline Facility bending-magnet beamline (BL20B) at the Photon Factory (PF), and for a Cu foil standard at the XAFS wiggler beamline at the Australian Synchrotron (AS). All XAFS data were recorded at 300 K, utilizing a conventional ion chamber configuration, and with harmonic rejection accomplished by detuning the monochromator. The χ(k) extraction and fitting procedures were performed using the Ifeffit package, with Artemis version 0.8.014, and Athena version 0.8.061 (Ravel & Newville, 2005). The Fourier transform operations required for the χ(k) to χ(R) conversion were performed using a Hanning window function (with Δk = 1.0 Å⁻¹). Further data transformation details are given in the footnotes to Tables 4 to 8.

3. Computational methods

3.1. MD simulations

MD simulations of metallic crystallites at 300 K were performed with the Kalypso package, version 3.1 (Karolewski, 2005). The arrangements and orientations of the atomic layers used to construct the simulation crystallites are summarized in Table 1. The lattice constants for Cu (Kroeger & Swenson, 1977), Ni (Collins & Gehlen, 1971), Fe (Basinski et al., 1955), Ti (Collins & Gehlen, 1971) and Au (Martienssen, 2005) refer to 300 K. Periodic boundary conditions were applied to the crystallite along the x and y directions (side faces). The crystallites were terminated by free surfaces in the z direction (top and bottom faces). Crystallite sizes much larger than the information range of XAFS (∼10 Å, determined by electron mean free paths of 20 Å or less) were employed, both to avoid coarse-sampling of reciprocal space, which could lead to an underestimation of thermal vibrational amplitudes (Winkler & Dove, 1992), and to avoid surface relaxation effects in the configurational data. Temperatures were maintained by means of a Berendsen thermostat (Berendsen et al., 1984). MD simulations were carried out using a 1 fs time-step. After equilibration for 10 ps, atomic configurational data were stored at 5 ps intervals until termination at 60 ps.

Table 1
Structures of metal crystallites used for MD simulations at 300 K

Metal	Number of layers	Atoms per layer	Layer orientation	Lattice constants (Å)
Cu (f.c.c.)	30	625	{100}	3.6149
Ni (f.c.c.)	30	625	{100}	3.5241
Fe (b.c.c.)	30	625	{100}	2.8665
Ti (h.c.p.)	60	693		a: 2.9505; c: 4.6830
Au (f.c.c.)	30	625	{100}	4.0784

MD simulations were performed using six EAM (or EAM-like) potentials for Cu, Ni, Fe and Au, and four EAM-like potentials for Ti. Table 2 lists the abbreviations used in this paper for each potential, and indicates the functional type of each potential. Most of the potentials have analytic forms. The TB potentials are distinguished in terms of whether they extend to the second (TB2) or fifth (TB5) neighbour distance. The Voter (Cu, Ni, Au), GRS (Au) and ZM (Ti) potentials were implemented using numerical tables provided by their originators.

Table 2
List of interatomic potentials for Cu, Ni, Fe, Ti and Au used in this study

The table refers to the potentials using the abbreviations used in the text. The potential functional types are also indicated.

Potential	Type†	Reference	Potential	Type†	Reference
Cu			Fe
TB5	TB	Cleri & Rosato (1993)	DD	EAM	Dudarev & Derlet (2007)‡
TB2	TB	Karolewski (2001)	OLS	EAM	Olsson (2009)
RTS	FS	Rafii-Tabar & Sutton (1991)	FSA	FS	Ackland et al. (2004)
Voter	EAM	Voter (1998)	EFS	EFS	Dai et al. (2006)
GTL	EAM	Gong et al. (2004)	FSM	FS	Marchese et al. (1988)
AV	FS	Ackland & Vitek (1990)	FSL	FS	Lau et al. (2007)
Ni			Au
TB5	TB	Cleri & Rosato (1993)	TB5	TB	Cleri & Rosato (1993)
TB2	TB	Karolewski (2001)	TB2	TB	Karolewski (2001)
ATVF	FS	Ackland et al. (1987)	AV	FS	Ackland & Vitek (1990)
EFS	EFS	Dai et al. (2006)	EFS	EFS	Dai et al. (2006)
VC	EAM	Voter & Chen (1987)	Voter	EAM	Voter (1993)
ZLL	EAM	Zhang et al. (1998)	GRS	EAM	Grochola et al. (2005)
Ti			Ti
ZM	EAM	Zope & Mishin (2003)	FSA	FS	Ackland (1992)
IKV	FS	Igarishi et al. (1991)	LREP	LREP	Dai et al. (2009)

†TB: tight-binding; FS: Finnis-Sinclair; EAM: embedded atom method; EFS: extended Finnis-Sinclair; LREP: long-range empirical potential.
‡Note that coefficients f₆ and V₅ are reported with incorrect signs therein.

Minor adjustments were made to some potentials. An interpolation function was used to truncate all abruptly terminating potentials smoothly in the region beyond the potential cut-off distance. The length scale parameter for the original RTS potential (Cu) was adjusted from 3.6100 to 3.6149 Å in order to reproduce the lattice constant of Cu at 300 K. The length scale of the OLS (Fe) potential was adjusted by a factor of 0.99875 in order to reproduce the lattice constant of Fe at 300 K (2.8665 Å). Other potentials were used without adjustment.

3.2. XAFS prediction and fitting methods

3.2.1. Cumulant predictions

Cumulants of XAFS scattering path distributions were extracted from the MD configurational data using a nested neighbour search procedure. All atoms selected for use in cumulant predictions and XAFS simulations (§3.2.2) were located in bulk-like sites that were situated at more than twice the potential cut-off distance from the free and periodic surfaces of the simulation crystallites. The low-order cumulants (C_n) and moments (μ_n) of the scattering path distribution, ρ(r), are related through the following expressions, which depend on the variable r (the instantaneous half-length of the scattering path, with mean value R),

Thus, R represents the mean scattering path half-length, while σ² represents the mean square variation of the scattering path half-length (Fornasini, 2001). For single-scattering paths, the MD cumulants reflect the distribution of the instantaneous relative displacements of two atomic centres. For single-scattering paths, R is the mean distance between the absorber and scattering atoms, while σ² is their mean square relative displacement. The first cumulant is conventionally described using a related quantity, ΔR,

where R₀ represents the scattering path half-length when the absorber and scattering atoms are located at their ideal lattice site positions. In this work, R₀ is always identical to the experimental crystallographic bond length (for single-scattering paths), or the sum of such bond lengths (for multiple-scattering paths). MD predictions of ΔR, σ² and C₃ will be presented later in this paper. Scattering paths are distinguished by the notation used by Binsted et al. (2005). For example, `0-2-0' designates a (two-leg) single-scattering path that involves the absorber atom (`0') and an atom in the second coordination shell (`2').

The cumulants arising from XAFS fitting procedures have a different significance from those predicted by MD. Whereas the MD cumulants represent the true moments of the scattering path distribution (the `real' distribution), ρ(r), the XAFS cumulants are the moments of a weighted or `effective' distribution, P(r),

that arises from the XAFS equation (λ is the photoelectron mean free path) (Bunker, 1983; Fornasini et al., 2004). The distinction between the real and effective distributions is mainly of importance for predictions of the first cumulant.

3.2.2. XAFS simulations

Synthetic XAFS spectra, i.e. the XAFS modulation functions χ(k), were calculated from the MD configurational data using the ab initio FEFF code, version 8.2 (Rehr & Albers, 2000). The calculation procedure for χ(k) entails the production and execution of a FEFF input file for each absorber site and simulation time. The resulting χ(k) includes contributions from all scattering paths with one to four legs up to a maximum half path length R_max = 6 Å. The total number of absorber atoms sampled for the χ(k) calculations ranged from 35000 to 80000 for the various metals. The averaged χ(k) for simulation times of 10–30 ps did not differ significantly from those for times of 35–60 ps. The site- and time-averaged χ(k) functions were adjusted to take account of the energy origin shift (E₀) and XAFS amplitude reduction factor (S₀ ²) deduced from experimental data, as discussed in §4.2. A normalization (McMaster) correction (Rehr et al., 1991) was also applied to χ(k). The conversion of the adjusted χ(k) function to the synthetic XAFS spectrum in R-space, χ(R), was performed with Athena (Ravel & Newville, 2005), using a Hanning window function (with Δk = 1.0 Å⁻¹).

4. Results

4.1. MD cumulant predictions

MD simulations can support two distinct approaches to the analysis of XAFS data: the MD cumulants can be extracted for use in an XAFS cumulant analysis procedure (Edwards et al., 1997; Bunker, 2010), or synthetic XAFS spectra can be predicted from MD configurational data (Binsted et al., 2005). In this section the MD cumulants predicted by different potentials for each metal are first compared.

Figs. 1 to 5 plot the first and second MD cumulants (ΔR and σ², respectively) predicted for single- and multiple-scattering paths of Cu, Ni, Fe, Ti and Au, that have scattering path half-lengths (R₀) up to 5 Å. The third MD cumulants (C₃), predicted for the 0-1-0 scattering paths only, are listed in Table 3, while Fig. 6 displays C₃ predictions for a range of scattering paths in Cu and Au.

Table 3
Third cumulants (C₃) predicted by MD models for the 0-1-0 paths of Cu, Ni, Fe, Ti and Au (units: 10⁻⁴ Å⁻³)

Experimental values obtained from the CD fitting models are also shown for Cu, Ni and Au.†

Cu	Model	TB5	TB2	GTL	RTS	AV	Voter	CD
	C3‡	1.5	2.1	1.8	4.1	0.7	1.5	2.2 ± 0.3
Ni	Model	TB5	TB2	ATVF	EFS	VC	ZLL	CD
	C3	0.9	1.4	0.7	0.8	0.7	0.1	0.5 ± 0.4
Fe	Model	DD	OLS	FSA	EFS	FSM	FSL
	C3	0.2	0.5	0.8	0.8	0.2	0.2
Ti	Model	ZM	IKV	FSA	LREP
	C3	1.0	1.1	1.0	1.7
Au	Model	TB5	TB2	AV	EFS	Voter	GRS	CD
	C3	4.7	3.7	2.5	4.4	4.3	5.1	1.2 ± 0.6

†Maximum statistical errors (10⁻⁴ Å⁻³) in the MD predictions of C₃ are Cu: ±0.09; Ni: ±0.03; Fe: ±0.04; Ti: ±0.06; Au: ±0.09.
‡Mean values from fits to APS, AS and PF data.

Figure 1 First and second cumulants (ΔR and σ^2, respectively) for several XAFS scattering paths predicted by MD simulations for Cu at 300 K. Data for multiple-scattering paths are indicated by arrows.

Figure 2 First and second cumulants (ΔR and σ^2, respectively) for several XAFS scattering paths predicted by MD simulations for Ni at 300 K. Data for multiple-scattering paths are indicated by arrows.

Figure 3 First and second cumulants (ΔR and σ^2, respectively) for several XAFS scattering paths predicted by MD simulations for Fe at 300 K. Data for multiple-scattering paths are indicated by arrows.

Figure 4 First and second cumulants (ΔR and σ^2, respectively) for several XAFS scattering paths predicted by MD simulations for Ti at 300 K. Data for multiple-scattering paths are indicated by arrows.

Figure 5 First and second cumulants (ΔR and σ^2, respectively) for several XAFS scattering paths predicted by MD simulations for Au at 300 K. Data for multiple-scattering paths are indicated by arrows.

Figure 6 Third cumulants (C3) for several XAFS scattering paths predicted by MD simulations for Cu (AS data) and Au at 300 K. Data for multiple-scattering paths are indicated by arrows.

Among different potentials, the predictions of ΔR for the 0-1-0 scattering paths (i.e. the paths with lowest R₀ in Figs. 1 to 5) are dispersed over a typical range of 0.01 Å or less. However, the ΔR predictions for the 0-1-0 paths based on the RTS potential for Cu, and the ZLL potential for Ni, deviate significantly from the respective group mean values. The ΔR predictions for different potentials tend to maintain their relative size order, but diverge from each other as R₀ increases.

With a few exceptions, the MD predictions of σ² for the 0-1-0 scattering paths show a tendency to cluster in a range ≤0.002 Ų. For longer paths, the range of predicted values tends to increase. The σ² estimates for some potentials (Cu: RTS; Fe: EFS; Ti: IKV; Au: AV) consistently deviate from their respective group mean values. The σ² predictions for multiple-scattering paths are usually close to those for the single-scattering paths that have comparable path lengths. The σ² predictions for different potentials also tend to maintain their relative size order as R₀ increases. However, this order is unrelated to the size order observed for the ΔR predictions.

The C₃ predictions produced by different potentials display only order-of-magnitude consistency (Fig. 6; C₃ predictions for Ni, Fe and Ti may be found in the supplementary material for this paper, Figs. S1 to S3¹). This behaviour is discussed in §5.1. For most metal/potential combinations, C₃ has its highest value for the first shell (Cu, Ni and Au) or for the first and second shells (Fe and Ti). All potentials for Cu, Ni and Au predict minimum C₃ values for the second shell (0-2-0 path). The ZLL potential for Ni predicts negative C₃ for all paths, whereas the predictions for other Ni potentials are typically positive for paths other than 0-2-0.

4.2. XAFS cumulant fitting

The MD cumulants predicted by different potentials were evaluated for use in XAFS analysis by employing them to fit experimental XAFS data for metal foil samples in a cumulant expansion model (Bunker, 1983, 2010). In the fitting procedure, the first (R) and second (σ²) MD cumulants of each path were optimally adjusted using global (path-independent) scale factors α and β, respectively. The scale factors were included in order to compensate for (i) intrinsic errors in the predicted MD cumulants, (ii) the systematic differences that are known to exist between the MD and XAFS cumulants (§3.2.1) and (iii) experimental measurement errors. Specifically, the first and second cumulants for the nth scattering path were fitted by making the following adjustments to the predicted MD cumulants: and (where R_n and are the first and second MD cumulants). In this scheme, = . For the 0-1-0 paths of the f.c.c. metals (Cu, Ni, Au), the third MD cumulant (C₃) was also included in the fit, with a corresponding scale factor γ.

Since the MD cumulants are scaled during the fitting procedure, the fitting error metric (R-factor) characterizes the ability of the MD simulations to predict the relative values of the optimum XAFS cumulants for a range of scattering paths. The fitted values of α, β and γ represent the correction factors that must be applied to the MD cumulants in order to produce the best fit with the experimental XAFS data.

An accurate potential would be expected to produce a small R-factor, with physically reasonable fitted values of α and β. Since = R_n-R₀, the fitted value of α is expected to be small but non-zero. The fitted value of β represents a scaling factor applied to the predicted second cumulants, and adjusts the latter to optimize the fit with the experimental second cumulant (if β = 1.0, no adjustment is required). It is worth emphasizing that the values of both α and β are also influenced by any measurement errors in the spectra to which they are fitted.

The values of E₀ and S₀ ² used for the MD cumulant fits were first estimated independently of the MD predictions by fitting the experimental XAFS data (single- and multiple-scattering paths) with a simple model. The first cumulant was fitted as = , where α is treated as a path-independent coefficient. The correlated Debye (CD) model (Sevillano et al., 1979) was used with a fitted Debye temperature, to estimate the second cumulants. The third cumulant was also fitted, for the first shell only. The estimated values of E₀ and S₀ ² were not found to be sensitive to the fitting scheme, and several different approaches yielded similar estimates. These fits intentionally covered a limited range of scattering paths (details in Tables 4 to 8). Since R₀ ≃ R (to within <0.1%), the α parameters in the MD and CD fitting models have a similar significance, although they are not strictly identical.

Table 4
R-factor and scaling parameters (α, β, γ) required for MD cumulants (R, σ², C₃ respectively), as obtained from fits of different MD models to XAFS data for Cu (300 K) measured at the AS, APS and PF†

Fitting results are also shown for the correlated Debye (CD) model.

		Model‡
Parameter	Data	TB5	TB2	GTL	RTS	AV	Voter	CD§
R-factor	AS	0.0030	0.0032	0.0030	0.0024	0.0032	0.0031	0.0047
	APS	0.0026	0.0029	0.0026	0.0021	0.0027	0.0027	0.0041
	PF	0.0050	0.0042	0.0045	0.0029	0.0050	0.0042	0.0077
α (× 10^4)	AS	3	3	3	−3	4	3	5
	APS	13	13	13	7	14	13	15
	PF	−45	−45	−45	−52	−44	−46	−44
β	AS	1.099	0.947	0.970	0.723	1.082	1.104	–
	APS	1.164	1.004	1.028	0.767	1.147	1.170
	PF	1.134	0.978	1.002	0.747	1.118	1.141
γ	AS	1.3	0.9	1.1	0.4	2.8	1.3	–
	APS	1.5	1.1	1.3	0.5	3.3	1.6
	PF	1.0	0.7	0.9	0.3	2.3	1.1

†Fitted in R-space, with k² weighting, over ranges: k = 3–16 Å⁻¹ (AS) or 3–15 Å⁻¹ (APS and PF), and R = 1.5–4.8 Å. All fits above use values of E₀ and S₀ ² that were estimated from prior fits to a CD model in the range R = 1.5–3.6 Å (AS data: E₀ = 5.4 eV, S₀ ² = 0.89; APS data: E₀ = 6.00 eV, S₀ ² = 0.92; PF data: E₀ = 6.60 eV, S₀ ² = 0.85). All paths with ≤4 legs and R_eff < 5.3 Å (16 paths in total) are included in the fits.
‡Typical fitting errors for the MD models are α: ±5 × 10⁻⁴; β: ±0.015; γ: ±0.2.
§Fitted Debye temperatures (±4 K) are 329 K (AS), 319 K (APS), 323 (PF).

Table 8
R-factor and scaling parameters (α, β, γ) required for MD cumulants (R, σ², C₃ respectively), as obtained from fits of different MD models to XAFS data for Au (300 K) measured at the APS†

Fitting results are also shown for the correlated Debye (CD) model.

	Model‡
Parameter	TB5	TB2	AV	EFS	Voter	GRS	CD§
R-factor	0.0111	0.0111	0.0113	0.0117	0.0112	0.0108	0.0579
α (× 10^4)	−39	−39	−39	−38	−39	−39	−37
β	0.601	0.680	0.779	0.644	0.606	0.568	–
γ	0.30	0.39	0.56	0.33	0.34	0.28	–

†Fitted in R-space, with k² weighting, over ranges: k = 2.0–15.5 Å⁻¹, and R = 1.8–5.2 Å. All fits above use values of E₀ and S₀ ² that were estimated from prior fits to a CD model in the range R = 1.8–3.5 Å (E₀ = 5.2 eV, S₀ ² = 0.90). All paths with ≤4 legs and R_eff < 5.8 Å (15 paths in total) are included in the fits.
‡Typical fitting errors for the MD models are α: ±7 × 10⁻⁴; β: ±0.01; γ: ±0.08.
§Fitted Debye temperature is 193 ± 3 K.

Tables 4 to 8 summarize the results of fitting the experimental XAFS spectra for Cu, Ni, Fe, Ti and Au with the MD cumulant models. The fitting procedures include all scattering paths that influence the XAFS spectra within the fitting ranges of ∼1.5 to 5 Å (details in Tables 4 to Table 7). Tables 4, 5, 6 and 8 (for the f.c.c. and b.c.c. metals) also include comparisons with a CD model that fits the XAFS spectra over the same range (and includes the third cumulant in the fits for the 0-1-0 paths of the f.c.c. metals). With a few exceptions (Ni: ZLL potential; Ti: IKV potential), the fitting capabilities of the optimally scaled MD cumulants for different potentials are found to be comparable for a given metal. The MD and CD models employ an equal number of fitting parameters (three for Cu, Ni, Au; two for Fe, Ti). On the basis of the R-factor, the scaled MD cumulant models typically fit the experimental XAFS spectra as well as, or better than, the CD models. Fig. 7 displays the results of fitting the MD cumulants predicted by different EAM potentials for Cu to the experimental χ(R) data, using the k-space weighting employed in the fits (details in Table 4). Similar graphs showing the fitting results for Ni, Fe, Ti and Au may be found in the supplementary material to this paper (Figs. S4 to S7).

Table 5
R-factor and scaling parameters (α, β, γ) required for MD cumulants (R, σ², C₃ respectively), as obtained from fits of different MD models to XAFS data for Ni (300 K) measured at the APS†

Fitting results are also shown for the correlated Debye (CD) model.

	Model‡
Parameter	TB5	TB2	ATVF	EFS	VC	ZLL	CD§
R-factor	0.0038	0.0044	0.0052	0.0040	0.0040	0.0108	0.0052
α (× 10^4)	1	1	2	1	1	4	2
β	1.214	0.939	1.154	1.041	1.177	1.068	–
γ	1.4	0.90	2.0	1.6	1.7	−12.7	–

†Fitted in R-space, with k² weighting, over ranges: k = 3–15 Å⁻¹ (APS), and R = 1.5–4.6 Å. All fits above use values of E₀ and S₀ ² that were estimated from prior fits to a CD model in the range R = 1.5–3.6 Å (E₀ = 10.1 eV, S₀ ² = 0.81). All paths with ≤4 legs and R_eff < 5.6 Å (16 paths in total) are included in the fits.
‡Typical fitting errors for MD models are α: ±5 × 10⁻⁴; β: ±0.025; γ: ±0.3 (for the ZLL model α: ±9 × 10⁻⁴; β: ±0.032; γ: ±4).
§Fitted Debye temperature is 407 ± 5 K.

Table 6
R-factor and scaling parameters (α, β) required for MD cumulants (R, σ², respectively), as obtained from fits of different MD models to XAFS data for Fe (300 K) measured at the APS.†

Fitting results are also shown for the correlated Debye (CD) model.

	Model‡
Parameter	DD	OLS	FSA	EFS	FSM	FSL	CD§
R-factor	0.0141	0.0153	0.0114	0.0195	0.0155	0.0156	0.0158
α (× 10^4)	−29	−30	−36	−33	−30	−30	−28
β	1.132	1.080	0.962	0.845	1.098	1.104	–

†Fitted in R-space, with k² weighting, over ranges: k = 3–15 Å⁻¹ (APS), and R = 1.5–5.0 Å. All fits above use values of E₀ and S₀ ² that were estimated from prior fits to a CD model in the range R = 1.5–4.0 Å (E₀ = 7.2 eV, S₀ ² = 0.92). All paths with ≤4 legs and R_eff < 5.7 Å (26 paths in total) are included in the fits.
‡Typical fitting errors for the MD models are α: ±5 × 10⁻⁴; β: ±0.04.
§Fitted Debye temperature is 426 ± 8 K.

Table 7
R-factor and scaling parameters (α, β) required for MD cumulants (R, σ², respectively), as obtained from fits of different MD models to XAFS data for Ti (300 K) measured at the PF†

	Model‡
Parameter	ZM	IKV	FSA	LREP
R-factor	0.0392	0.0427	0.0334	0.0331
α (× 10^4)	−57	−52	−54	−48
β	1.018	1.301	0.918	0.978

†Fitted in R-space, with k³ weighting, over ranges: k = 2–14 Å⁻¹, and R = 1.9–5.2 Å. All fits above use values of E₀ and S₀ ² that were estimated from prior fits to a CD model in the range R = 1.9–3.6 Å (E₀ = 7.4 eV, S₀ ² = 0.68). All paths with ≤4 legs and R_eff < 5.8 Å (24 paths in total) are included in the fits.
‡Typical fitting errors for MD models are α: ±1 × 10⁻³; β: ±0.05.

Figure 7 Symbols: experimental XAFS spectrum for Cu (measured at the AS); solid lines: XAFS spectra fitted using scaled MD cumulants (see text for fitting details).

In Table 4, the variation of α between different Cu potentials (6 × 10⁻⁴) is an order of magnitude smaller than the variation of α between different beamlines (4.8 × 10⁻³). Both the sign and the magnitude of α vary with synchrotron source, e.g. in the range −0.0044 to 0.0004 for the CD fitting model. A reviewer of this paper has pointed out that discrepancies of this kind might arise from various angle-to-energy conversion errors in the XAFS measurements that are not necessarily removable by shifting ΔE. The range of values obtained in this work for α is not unusual for recent studies of Cu (at 300 K) reported by other groups. For example, a study of Cu XAFS spectra from 11 synchrotron sources fitted α in the range −0.002 to 0.000 (three-shell fits) (Kelly et al., 2009). First-shell fits to Cu XAFS spectra by Newville et al. (2009) fitted α between −0.0073 and −0.0050 for a single beamline using different fitting procedures. Four-shell fits to Cu XAFS spectra from three established beamlines fitted α between −0.002 and −0.0004 (Gaur et al., 2013). The elimination of beamline-dependent errors in XAFS remains an active field of investigation (Chantler et al., 2012). In this study the variation in α between different beamlines contributes an instrumental error of about ±0.005 Å to the estimate of the Cu first shell distance.

For metals other than Au, the fitted values of β (the scale factor for the second MD cumulant) are scattered around a mean value near 1.0. For Au, the fitted values of β are always <1.0, with a mean value of 0.65. In Table 4, the variation of β among different Cu potentials (±19% of the mean value) is much larger than the variation of β between different beamlines (±3% of the mean value). The relative values of β predicted by different Cu potentials are beamline-independent (e.g. the TB2:TB5 ratio of β values is 0.8622 ± 0.0005 for all beamlines). This suggests that the small variation in β between different beamlines is due to experimental error rather than structural differences (e.g. static disorder) in the unannealed foil standards.

The fitted values of γ (the scale factor fitted for the third MD cumulant of Cu, Ni and Au only) span a relatively broad range that reflects the scatter in the corresponding MD cumulants (§4.1). The large relative errors in the estimates of γ (typically ±20%) are due to the high correlation (r > 0.7) that exists between γ and α. In Table 4, the variation of γ among different Cu potentials (±65% of the mean value) is somewhat larger than the variation of γ between different beamlines (±20% of the mean value). To a good approximation, the relative values of γ predicted by different Cu potentials are beamline-independent (e.g. the TB2:TB5 ratio of γ values is 0.71 ± 0.02 for all beamlines). This suggests that the scatter in the γ values between different Cu potentials is a real effect, rather than an artefact of the fitting procedure. For Ni with the ZLL potential, α and γ are anti-correlated (r = −0.76), in contrast to other f.c.c. models. For this potential, both fitted parameters display anomalous values that are suggestive either of a systematic error in the pair correlation function predicted by the ZLL potential, or fitting errors caused by interactions between α and γ in the fitting model.

For the metals other than Cu, the fitted values of α, β and γ found in Tables 5–8 refer to a single beamline and thus provide no information about any beamline-dependent error in the experimental data. To address this, Table 9 compares R and σ² values for one (Cu, Ni, Au) or two (Fe, Ti) coordination shells obtained in this work using the CD model whose fitted data are summarized in Tables 4 to 8, with representative fitted data derived from recent studies (only limited data for Ti are available at 300 K). Although the CD model was constrained to use only two (b.c.c., h.c.p.) or three (f.c.c.) fitting parameters over multiple coordination shells, its results are in reasonable agreement with the previous studies. The maximum discrepancy in R of 0.015 Å is observed for the Au data, while the σ² fits agree to within the fitting errors with one exception (Fe first shell).

Table 9
Comparison of fitted first and second cumulants (R and σ²) for one (Cu, Ni, Au) or two (Fe, Ti) coordination shells obtained in this work using the correlated Debye model (see Tables 4 to 8 for fitting procedures), with representative fits reported in the recent literature (Cu: Kelly et al., 2009; Ni: Krayzman et al., 2009; Fe: Gordon & Crozier, 2006; Ti: Felderhoff et al., 2004; Au: Comaschi et al., 2009)

	Present work		Literature
Metal	R (Å)	σ^2 (10^−4 Å^2)	R (Å)	σ^2 (10^−4 Å^2)
Cu	2.557 ± 0.002 (AS)	89 ± 2 (AS)	2.553 ± 0.003	95 ± 2
	2.560 ± 0.002 (APS)	93 ± 2 (APS)
	2.545 ± 0.002 (PF)	91 ± 3 (PF)
Ni	2.493 ± 0.002	65 ± 2	2.485 ± 0.002	64 ± 2
Fe	2.475 ± 0.002	63 ± 2	2.474 ± 0.006	52 ± 4
	2.858 ± 0.002	70 ± 3	2.852 ± 0.009	67 ± 8
Ti	2.878 ± 0.004	83 ± 5	2.920 ± 0.005†
	2.932 ± 0.004	85 ± 5
Au	2.870 ± 0.004	83 ± 4	2.8849 ± 0.0007	83.4 ± 0.3

†Average of first two shells.

4.3. Synthetic XAFS spectra

The relative accuracy of synthetic XAFS spectra can be expected to correlate, to a large extent, with the accuracy of the corresponding MD cumulant predictions. The FEFF calculations automatically assign the correct weights to each scattering path contribution, so the distinction between the real and effective distributions [equation (5)] does not have to be considered when evaluating synthetic XAFS spectra. The first cumulants (ΔR in Figs. 1 to 5) determine the positions of peaks in χ(R). The distribution of ΔR values for any given metal is relatively narrow for the majority of potentials (in a range ±0.005 Å for most 0-1-0 paths). The MD cumulant predictions are more clearly differentiated in terms of their ability to predict the second cumulants. The fitting parameter β (Tables 4 to 8) optimally scales the MD second cumulants (σ²) that determine the peak widths in χ(R). Values of β near 1.0 imply a greater similarity between the predicted MD second cumulants and the XAFS second cumulants. If β > 1, the predicted second cumulants will be too small, so the peaks in χ(R) will be too narrow, and vice versa. For the Au potentials, the fitted values of β range from 0.57 to 0.78 (Table 8), and thus the peak widths in synthetic XAFS spectra will be too high for all Au potentials. For the other metals (Cu, Ni, Fe and Ti) there is at least one potential for which β lies in the range 1.00 ± 0.05.

The most realistic synthetic χ(R) functions for Cu, Ni, Fe and Ti (as identified on the basis of the fitted β parameters) are compared with experimental XAFS χ(R) data in Fig. 8 (for Cu, the AS experimental data are selected). The synthetic χ(R) functions employ the same values of E₀ and S₀ ² as the fits to the cumulant expansion models (details in Tables 4 to 8). The unadjusted synthetic χ(R) are inferior in quality to the fits provided by the cumulant expansion model based on scaled MD cumulants (§4.2). However, for Fe (FSA potential), the synthetic and experimental χ(R) functions are in fair agreement up to 5 Å. For Cu (GTL potential), Ni (EFS potential) and Ti (LREP potential), the synthetic χ(R) data tend to display inaccuracies in peak positions and/or shapes above 3 Å. Binsted et al. (2005) performed direct XAFS simulations for bulk Cu based on use of the TB5 potential, but found limited agreement of peak heights, except for the 0-1-0 path. Higginbotham et al. (2009) did not report χ(R) predictions, but observed that the form of EAM-type potentials could significantly modify the χ(k) predictions for Fe. Taking into account the beamline-dependence of the fitted β values for Cu, it is not straightforward to decide which Cu potential produces the most accurate predictions of peak widths in χ(R), as judged by the proximity of β to 1.0. The GTL potential is most accurate for the AS and PF XAFS data, while the TB2 potential is most accurate for the APS XAFS data.

Figure 8 Symbols: experimental XAFS spectra for Fe, Ti, Cu (AS data) and Ni; solid lines: XAFS spectra predicted from MD configurational data for the same metals, using the indicated potentials.

5. Discussion

5.1. Reliability of MD cumulant predictions

5.1.1. First cumulants

MD simulations predict cumulants of the real path length distribution (§3.2.1). The first MD cumulant, R, differs significantly from the site–site distance, R₀, due to the effects of atomic thermal motion. The instantaneous relative thermal displacements, = , of atoms at two sites (labelled 1 and 2) can be resolved into components and that are oriented parallel and perpendicular, respectively, to the site–site vector (Fornasini et al., 2004). The deviation of the first MD cumulant from R₀ is approximately related to as follows,

If required, and can be predicted by projecting the current internuclear vector on the site–site vector (Sanson, 2010) (however, angular momentum is not conserved under periodic boundary conditions, so the prediction of requires a simulation method that explicitly eliminates errors due to rotational drift).

Experimental first-shell values estimated for δ at 300 K are 0.0039 Å for Cu (Fornasini et al., 2004) and 0.0011 Å for Au (Comaschi et al., 2009). These values will increase at higher temperatures. An accurate MD simulation should provide an estimate of the first real cumulant such that R ≃ , and thus ≃ δ. However, the first MD cumulant predictions (Figs. 1 to 5) are typically larger, such that ΔR = 0.01–0.02 Å for many 0-1-0 paths (e.g. 0.008 to 0.024 Å for Cu, and 0.017 to 0.026 Å for Au). These larger shifts are artefacts that originate from the practice of fitting interatomic potentials to static lattice properties, without regard for the lattice expansion that takes place after thermalization (Sheng et al., 2011). Specifically, simulation crystallites that are constrained to be periodic in two dimensions will relieve thermal stress via an expansion in the third, unconstrained, dimension. As a result, the MD prediction of ΔR is artificially high for any scattering path that has legs with a vector component lying normal to the free surface. These small structural distortions, due to anisotropic thermal expansion, cannot easily be avoided if the surface–vacuum interface is to be retained in the MD simulation. Scattering paths that lie entirely in the periodically constrained plane do not show this dispersion, but such paths are relatively rare, e.g. the 0-4-0 path for Ti (R₀ = 4.683 Å in Fig. 4), which lies parallel to the c axis.

5.1.2. Second cumulants

The accuracy of predicted thermal displacements is necessarily related to the ability of MD to model both (i) phonon occupation statistics, and (ii) phonon dispersion behaviour, at the temperature of interest. These are distinct issues.

MD ensembles reflect the properties of particle systems that obey classical (Boltzmann) statistics. However, in real metals, phonon occupation numbers, and thus properties that depend on phonons, are determined by quantum (Bose–Einstein) statistics (Cahill et al., 2003; Turney et al., 2009). The classical and quantum occupation number schemes begin to converge around the Debye temperature, θ_D, which is indicative of the thermal energy required to excite the high-frequency modes in the phonon density of states. For Au, θ_D = 170 K, while for other metals examined in this study θ_D takes values of 343 K (Cu) and 420–467 K (Ni, Ti, Fe). On this basis, the lattice dynamics of Au can be described by a classical model at 300 K, whereas Cu is a borderline case. At 300 K, the phonon occupation numbers for the high-frequency modes of Ni, Ti and Fe are significantly influenced by quantum statistics, although it should be noted that their phonon contributions to the specific heat are all within 10% of the classical value (i.e. 3k_B). Below θ_D, MD simulations always underestimate the thermal displacements. At 0 K, classical particles become stationary at their lattice sites, whereas quantum particles retain their zero point energies and associated vibrational displacements.

For metals with cubic structures (e.g. Cu, Ni, Fe and Au), EAM potentials can often reproduce the experimental elastic constants (C₁₁, C₁₂, C₄₄) with a precision of 1–2%. This implies that these potentials are able to reproduce phonon dispersion curves near the Brillouin-zone origin, because the elastic constants involve phonons near the zone origin. However, EAM potentials often show poor reproduction of the phonon dispersion curves near the Brillouin-zone boundaries. This is particularly true for potentials with a small number of fitting parameters (Bian et al., 2008). Phonon dispersion curves have been predicted for only a few of the potentials used in this study. The TB5 potentials overestimate the phonon cut-off frequencies for Cu and Ni by 5% and 16%, respectively, but underestimate that for Au by 30% (Cleri & Rosato, 1993). This has the effect of narrowing the phonon density of states in Au, leading to excessive thermal displacements (Kallinteris et al., 1997), and is a probable explanation for the consistently low values of β fitted to the MD second cumulants for Au in Table 8. The tendency for EAM-type potentials to overestimate thermal displacements in MD simulations of bulk Au (Kallinteris et al., 1997; Chamati & Papanicolaou, 2004) and Au nanoparticles (Roscioni et al., 2011) has been noted previously. The underlying cause is probably the influence of non-central many-body forces in Au that cannot be modelled using simple EAM functional forms (Cleri & Rosato, 1993; Bian et al., 2008). The phonon cut-off frequency for Cu is also substantially (20%) underestimated by the RTS potential (Kimura et al., 1998), which accounts for the unusually low value of β fitted for this potential.

The choice of potential cut-off distance in the fitting procedure has a minor influence on the predicted bulk material properties, but significantly affects the forces that act between neighbouring atoms, and hence the corresponding vibrational dynamics. Cut-off distances are arbitrarily employed in MD simulations for computational efficiency, and there are presently no clear criteria for fitting them optimally (Baskes et al., 2001). Similar remarks apply to the interpolation functions that are necessary to truncate some analytic potentials smoothly at the cut-off distance. The TB5 and TB2 potentials (Cleri & Rosato, 1993; Karolewski, 2001) used in this work for Cu, Ni and Au are fitted to similar material properties, and have identical (exponential) functional forms, but differ in terms of their cut-off distances, which lie above the second and fifth nearest neighbour distances, respectively. This difference is sufficient to produce a 10–30% variation in the predicted second cumulants.

Thermal lattice distortions also influence the MD predictions of σ², to the extent of increasing mean values by typically 5% (this can be established from detailed analysis of the directional behaviour of the σ² predictions). For the first shell of Cu, the MD σ² predictions obtained in this work using the TB5 potential differ by <1% from those obtained for bulk Cu in NVT Monte Carlo simulations by a Beccara et al. (2003), but are 7% higher than those reported by Binsted et al. (2005) in NVE MD simulations. Fitting errors in the experimental determination of σ² by XAFS for the first Cu shell are typically 5% (Newville et al., 2009).

5.1.3. Third cumulants

For most metals at 300 K, the third cumulants are small and XAFS data can be fitted reasonably well without using them. This is fortunate, because the prediction of third cumulants represents a severe test for MD simulations. The values of the third cumulants are associated with anharmonic thermal vibrations in lattices, which ultimately determine the third-order elastic constants. However, none of the potentials used in this study was explicitly fitted to third-order elastic constants. Chantasiriwan & Milstein (1996) argue that an EAM-type potential should include oscillatory terms, and be fitted to at least three coordination shells, for accurate modelling of third-order elastic constants. Few currently used EAM-type potentials meet these criteria. Thermal expansion artefacts (see above) might also be expected to unphysically skew the distribution of atomic positions, and thus artificially inflate any predicted values of C₃. However, no evidence was found for this. The C₃ predictions obtained with the TB5 potential for the first eight scattering paths of Cu (Fig. 6) are in close agreement with those reported by Binsted et al. (2005) for bulk Cu using the same potential. Some reasonable predictions of C₃ (in the sense that γ ≃ 1) do arise, probably fortuitously, from several of the potentials used for Cu, Ni and Au (Tables 4, 5 and 8).

5.2. Prospects for XAFS analysis

After the energy origin shift (E₀) and XAFS amplitude reduction factors (S₀ ²) have been established for the sample of interest (by fitting XAFS data for experimental standards), it is relatively straightforward to incorporate MD results in XAFS data analysis, either via a cumulant expansion model (§4.2) or by direct simulation of XAFS spectra (§4.3).

In this study, global scaling parameters α, β and γ have been used to incorporate the first, second and third MD cumulants, respectively, into XAFS data fitting procedures based on standard XAFS data analysis software (Ravel & Newville, 2005) and a cumulant expansion model. The scaling parameters compensate for the inherent differences between the real and effective distributions, and for any errors that might exist in the MD predictions and in the experimental data. In Cu, the fitted value of α is largely determined by beamline-dependent factors (e.g. experimental errors), whereas the beamline influence on the fitted values of β and γ is somewhat weaker. Therefore, for Cu (and possibly for other metals), β and γ are mainly determined by the properties of the potential (i.e. MD prediction errors).

For systems that display high structural disorder (e.g. clusters, irradiated materials), a cumulant expansion fitting model is not suitable (Bunker, 1983, 2010). Under these circumstances, XAFS data can in principle be analysed by coupling a heuristic search procedure, such as reverse Monte Carlo simulation (McGreevy & Pusztai, 1988; Di Cicco & Trapananti, 2005) or genetic algorithms (Dimakis & Bunker, 2006), to a multiple-scattering code that generates the XAFS signal for trial structures. In contrast to the purely numerical optimization performed by heuristic techniques, MD simulations of χ(R) can provide a first-principles approach to XAFS analysis in cases where it is possible to specify (actual or hypothetical) initial conditions for the system of interest. Unlike predicted MD cumulants, the fitting properties of direct simulations of χ(R) are not influenced by differences between the real and effective distributions, but they remain similarly sensitive to MD simulation errors and experimental errors. The references given in §1 provide examples of MD applications in XAFS analysis. In particular, Binsted et al. (2005), Roscioni et al. (2011) and Price et al. (2012) discuss fitting strategies (for nanoscale and bulk materials) in some detail.

6. Conclusions

This study evaluates the ability of MD simulations to support the analysis of XAFS data for metals. The low-order MD cumulants (ΔR, σ², C₃) for XAFS scattering paths were calculated for the metals Cu, Ni, Fe, Ti and Au at 300 K using 28 interatomic potentials of the EAM type. The MD cumulant predictions were evaluated within a cumulant expansion XAFS fitting model, using global (path-independent) scaling factors. Direct simulations of the corresponding XAFS spectra, χ(R), were also performed using MD configurational data in combination with the FEFF ab initio code (Rehr & Albers, 2000).

The scaling parameters that are fitted in the cumulant expansion model compensate for differences between the real and effective scattering path distributions, and for any errors that might exist in the MD predictions and in the experimental data. The fitted value of ΔR is particularly susceptible to errors that arise both from experimental factors and in­advertent lattice thermal expansion in the simulation crystallites. The unadjusted predictions of σ² vary in accuracy, but do not show a consistent bias for any metal except Au, for which all potentials overestimate σ². The unadjusted C₃ predictions produced by different potentials display only order-of-magnitude consistency. The suitability of current EAM-type potentials for predictions of C₃ is questionable, since their fitting databases do not include third-order elastic constants.

The accuracy of direct simulations of χ(R) for a given metal varies among the different potentials. For each of the metals Cu, Ni, Fe and Ti, at least one of the tested potentials was found to provide a reasonable simulation of χ(R), in the sense of reproducing the major peak positions and widths in the region R < 5 Å (Fig. 8). However, none of the potentials tested for Au was able to reproduce the peak widths with sufficient accuracy to be of value for XAFS data analysis.

It is difficult to anticipate how a given interatomic potential will perform when used for XAFS applications. MD predictions of ΔR could be improved if the lattice expansion due to thermal motion were included explicitly in the fitting procedure (Sheng et al., 2011). Accurate predictions of σ² require potentials that can reproduce phonon dispersion properties (in addition to elastic constants). Even then, quantum statistical effects may limit the accuracy of σ² predictions below the Debye temperature. The potential cut-off distance may also be a significant fitting parameter for XAFS applications of MD.

The systematic differences between the cumulants of the real and effective scattering path distributions complicate the fitting of XAFS data by MD cumulants. However, it is possible to compute cumulants for the effective distribution. Such calculations require a knowledge of the electron mean free path dependence on the wavevector, λ(k), which can be obtained from FEFF or other theoretical XAFS codes. Theoretical relationships between the real and effective cumulants have been derived that can be used for consistency checks (Bunker, 1983, 2010).

Footnotes

¹Supplementary data for this paper are available from the IUCr electronic archives (Reference: HF5227 ). Services for accessing these data are described at the back of the journal.