Chapter 5.14. Reverse Monte Carlo and molecular-dynamics approaches to EXAFS analysis

The problem of quantitative extended X-ray absorption fine-structure (EXAFS) analysis has been attracting attention for many years (Stern, 1974; Brouder et al., 1989). A common approach relies on the multiple-scattering (MS) formalism (Rehr & Albers, 2000; Natoli et al., 2003), where EXAFS is described by an infinite series of scattering events in a many-atom system. The fast convergence of the MS series is ensured by the finite lifetime of the excitation as well as static and thermal disorder in the material. Nevertheless, the accurate evaluation of all required MS contributions with proper damping amplitudes due to disorder effects is still a challenging task. Conventional EXAFS data fitting to experimental or theoretical standards (Lee et al., 1981) is therefore limited to analysis of the nearest coordination shells around the absorbing atoms only, since (i) the number of structural parameters [interatomic distances, mean-square relative displacement (MSRD) factors etc.] that need to be included in the analysis increases exponentially with an increase in the number of atoms in the neighbourhood of the absorber and (ii) the contributions of distant scattering paths usually overlap and thus the estimators of structural parameters are strongly correlated (Kuzmin & Chaboy, 2014). Moreover, conventional EXAFS analysis requires the a priori assumption of some shape (often Gaussian) for the distribution of inter­atomic distances in the material. Hence, even for the first coordination shell this can be inaccurate for a significantly distorted environment or disordered systems, where such a distribution can often be far from the Gaussian distribution (Filipponi, 2001; Clausen & Nørskov, 2000; Chill et al., 2015). Therefore, new methods of EXAFS analysis are required to treat these issues.

A possible solution is to rely on progress in simulations of material structures and dynamics. It has been shown in recent decades that theoretical modelling is able to complement the information available from experiments and that their combined use can also help to gain additional information that is not easily accessible from analysis of the experimental data only (McGreevy, 2001; Billinge & Levin, 2007; Prasai et al., 2015).

The best-known simulation-based approaches employed for EXAFS data analysis are the reverse Monte Carlo (RMC) and molecular-dynamics (MD) methods.

As in many other mathematically incorrect problems, one needs to introduce some constraints to stabilize the model with many degrees of freedom that is required to parameterize EXAFS spectra. A natural set of constraints can be obtained by acknowledging that the sought structural parameters are not actually independent: the interatomic distances in the distant coordination shells are geometrically linked to the interatomic distances in the nearest coordination shells, the distributions of bonding angles are closely related to the distributions of atomic deviations etc. To benefit from this purely geometrical information, one can reformulate the problem so that instead of directly fitting the average inter­atomic distances, MSRD factors etc., one can use a set of atomic coordinates in a sufficiently large piece of the investigated material as unknown variables of the model. The goal of the RMC procedure is to find such a set of atomic coordinates so that the corresponding calculated properties of the material (for example its EXAFS spectrum) are as close as possible to the experimental data.

A scheme of the conventional RMC approach, as proposed by McGreevy and Pusztai in the 1980s (McGreevy & Pusztai, 1988; Keen & McGreevy, 1990), is given in Fig. 1. Originally, it was designed for the interpretation of X-ray and neutron scattering data, but almost immediately its potential for EXAFS data analysis was also realized (Gurman & McGreevy, 1990).

Figure 1 Flowchart of the reverse Monte Carlo method for EXAFS data analysis (Gurman & McGreevy, 1990).

RMC simulation starts with an arbitrary initial configuration of atoms in a cell with a chosen size and shape. Periodic boundary conditions can be applied if necessary to avoid surface-related effects. Most RMC implementations do not allow one to vary the number of atoms during the simulation; hence, the density of the material should be known in advance and the volume of the cell and the number of atoms should be chosen appropriately. Note that these rather obvious requirements are not always easily met for disordered systems (Walters & Newport, 1996). For the modelling of crystalline materials, a natural choice for the initial configuration is the equilibrium structure, which can be obtained from diffraction studies. Note that the structural parameters of a crystalline material can be determined by diffraction techniques with much higher accuracy (better than 10⁻³ Å) than that provided by modern EXAFS analysis (usually about 10⁻² Å). Therefore, one may want to keep the size of the simulation cell fixed during the RMC–EXAFS modelling of crystalline materials, thus indirectly accounting for the information available from diffraction data even without simulation of the diffraction pattern.

At each RMC iteration, a random modification of the constructed atomic configuration is performed. Depending on the implementation and the investigated problem, such modification can involve small random displacements of all or a randomly selected few atoms, the swapping of two randomly picked atoms etc.

For the new proposed atomic configuration, the corresponding EXAFS spectrum is then calculated using ab initio codes such as FEFF (Ankudinov et al., 1998; Kas et al., 2024) or GNXAS (Filipponi et al., 1995; Filipponi, Di Cicco et al., 2024; Filipponi, Natoli et al., 2024). Typically, there are several absorbing atoms in the modelled atomic configuration. Independent EXAFS calculations should be carried out for each of them and the obtained results should be averaged to obtain a single configuration-averaged EXAFS spectrum, which then can be compared with the experimental EXAFS data. Importantly, the RMC approach allows one to benefit from the information encoded in the multiple-scattering contributions to the total EXAFS spectra; hence, these contributions should also be included in the analysis (Di Cicco et al., 2014).

The difference between the experimental and simulated data is characterized by a single number ξ. It can be defined, for instance, as a squared Euclidean distance between the total calculated (χ^calc) and experimental (χ^exp) EXAFS in k-space or R-space (or simultaneously in k-space and R-space if the wavelet transform is used; Timoshenko et al., 2012). The proposed atomic configuration will then be either accepted or discarded depending on the value of ξ. It is important to emphasize that the ξ value can be calculated for a single experimental EXAFS spectrum as well as for several EXAFS spectra acquired at different absorption edges. Moreover, additional experimental information can be used. In this case , where ξ_j is the difference between experiment and calculation for the jth data set and w_j are the corresponding weights. Commonly, for instance, the agreement between simulated and experimental X-ray and neutron scattering data is used together with EXAFS data (Bridges et al., 2014). Other structure-related properties can also be employed in this step: Bragg diffraction profiles, electron diffraction patterns, electric polarization values, the intensities of some X-ray absorption near-edge structure (XANES) features (Levin et al., 2013), information from NMR spectra (Wicks et al., 1997), anomalous X-ray scattering data (Hosokawa et al., 2011) etc. The disagreements between various experimental and calculated properties can provide a detailed description of the goodness of the proposed atomic configuration. Also, various chemical and geometrical constraints can be implemented at this stage by assigning some penalty (an increase in the ξ value) when the atoms become too close or too far from each other, when nonphysical values of some specific bonding angle are encountered (Tucker et al., 2007), when the number of nearest neighbours for some atom deviates from that expected (McGreevy, 2001) etc. The ability of the RMC method to benefit from such diverse experimental and theoretical information is one of the particular strengths of this method. On the other hand, one should be aware that only in rare cases can one ensure that the sample and experimental conditions are actually the same during, for example, XAS and electron diffraction measurements. Experimental data provided by different techniques may have quite different noise levels and different systematic errors, and the accuracy of the theories used for the calculation of different material properties may also vary significantly. Therefore, the assignment of proper weights w_j is often a non­trivial task.

The decision whether to accept or discard the proposed atomic configuration is made based on the so-called Metropolis algorithm (Metropolis et al., 1953). The value of ξ calculated for the proposed atomic configuration is compared with that obtained for the previous atomic configuration and is further denoted ξ^old. If ξ < ξ^old then the proposed atomic configuration is better than the previous configuration and should be accepted. However, if one discards all atomic configurations for which ξ > ξ^old, the difference ξ will always decrease and after some number of steps it will reach the local minimum. In order to ensure that the configurational space is sampled properly and that the global minimum can be found, the Metropolis algorithm also accepts some atomic displacements for which ξ > ξ^old with a probability proportional to exp[(ξ^old − ξ)/Θ], where Θ is a scaling parameter. For sufficiently low values of the parameter Θ one may expect that after a large number of iterations the modelled atomic configuration will be close to that which gives the minimal possible value of ξ (Narayan & Young, 2001). The value of the scaling parameter Θ is often related to the level of experimental noise (McGreevy, 2001). Alternatively, a so-called simulated-annealing approach (Kirkpatrick et al., 1983) can be applied: in this case the parameter Θ is not fixed during the simulations but decreases slowly (Timoshenko et al., 2012).

Knowing the atomic coordinates for the final atomic configuration one can calculate a set of structural parameters of interest, such as radial distribution functions (RDFs), average interatomic distances, atomic displacements etc. Since the RMC method relies on a stochastic process, the final solution (the final set of atomic coordinates) will be different if calculations are repeated several times starting from different initial configurations or employing different sequences of pseudorandom numbers. The results of these calculations are nevertheless expected to be close from a statistical point of view: RDFs, bonding angles and MSRDs will be close provided that the experimental data used contain a sufficient amount of information to recover these structural parameters. If the information content in the used set of experimental data is not sufficient to reconstruct the structure unambiguously, the RMC method tends to converge to the most disordered solution of those that are consistent with the provided experimental data (Tucker et al., 2007).

One of the most significant limitations of the RMC approach is that it is computationally inefficient. Many thousands of iterations are typically required before a reasonable agreement between experiment and theory is obtained. Moreover, the required computational time increases exponentially with an increase in the size of the configurational space that the algorithm needs to explore (for example with an increase in the maximal allowed atomic displacement) and with an increase in the dimensionality of the problem (for example with an increase in the number of atoms). Therefore, alternative optimization schemes have been suggested. For instance, it has been shown that the implementation of a so-called evolutionary algorithm (Holland, 1992) in the conventional RMC approach can reduce the required computational time significantly, especially for EXAFS analysis in crystalline systems (Timoshenko, Kuzmin et al., 2014).

Some software packages for RMC–EXAFS simulations include RMCProfile (Tucker et al., 2007), RMC++ (Gereben et al., 2007) and EvAX (Timoshenko, Kuzmin et al., 2014).

Unlike the RMC method, the molecular-dynamics (MD) approach yields a time-dependent 3D model of the investigated structure. In both classical and ab initio MD approaches the initial structure model evolves in time according to classical Newtonian laws of motion, which allow one to account for the atomic thermal disorder in the system. Note that since both methods describe the motion of atomic nuclei classically, such MD approaches cannot be used to model the motion of atoms at low temperatures, where the zero-point oscillations of atoms play an important role. To treat this issue, methods such as path-integral Monte Carlo (a Beccara et al., 2003) or path-integral MD (Berne & Thirumalai, 1986) should be used instead. On the other hand, at high temperatures MD methods usually perform well, and can be employed, for instance, to model anharmonic motion of atoms, leading to asymmetric distributions of bond lengths.

In ab initio MD the Born–Oppenheimer approximation is commonly employed, i.e. it is assumed that the relatively slow motion of heavy atomic nuclei can be separated from the motion of electrons. Therefore, the Schrödinger equation is solved for electrons in the electrostatic field produced by a frozen configuration of atomic nuclei for the calculation of forces. Typically, such calculations are carried out within density-functional theory (DFT) formalism. The contributions of core electrons are usually defined using some pseudo­potential approximation, while the valence electrons are treated explicitly, using, for example, the plane-wave approach (Ferlat et al., 2005). Once the potential energy V(r₁, r₂, …, r_n) of the system is known, one can calculate the force F_i acting on the ith atom as , where the vectors r₁, r₂, …, r_n correspond to the positions of atomic nuclei. Since the forces need to be recalculated at each iteration, such an approach (sometimes called Born–Oppenheimer MD) is obviously extremely computationally expensive. In some cases the efficiency of ab initio MD can be improved using the so-called Car–Parrinello method (Car & Parrinello, 1985; Spezia et al., 2006).

In classical MD the potential energy function V(r₁, r₂, …, r_n) is defined empirically (Ferlat et al., 2005), thus significantly reducing the requirements for computational resources. The use of the classical MD approach hence allows one to model much larger systems for longer periods of time (Delaye, 2001). Systems with millions of atoms and on timescales of up to milliseconds can currently be simulated, and complex phenomena such as phase transitions, self-assembly of molecular systems, protein folding and various biomolecular functionalities can be observed in silico (Karplus & McCammon, 2002; Klein & Shinoda, 2008). The accuracy of the obtained results is determined by the accuracy of the empirical force-field model.

The design of an empirical force field is an art in itself. Its functional form and a set of parameters should be chosen to reproduce the structure of the material and its known properties during simulations. A force-field model should ensure the correct coordination of atoms, bond lengths and lattice constants. Often, the fitting of elastic and dielectric properties is used to find the values of force-field parameters (Gale & Rohl, 2003). Additionally, phonon spectra, known either from experiment or from ab initio calculations (Kuzmin & Evarestov, 2009), can be used to fit the values of force-field parameters. An alternative approach for force-field model optimization is fitting of the parameterized V(r₁, r₂, …, r_n) function to the potential energy function derived directly from ab initio simulations, thus allowing one to implicitly include information from DFT modelling in the classical MD simulations (Gale & Rohl, 2003).

The first attempts to correlate the information from MD simulations with the results of EXAFS analysis date back to the 1980s (Hayes & Boyce, 1980). Attempts to use MD for direct simulation of EXAFS spectra, however, are more recent (Palmer et al., 1996). The basic scheme of such MD–EXAFS analysis is given in Fig. 2. The MD models that are useful for EXAFS data analysis can be quite simple, since the MD approach is often just used to incorporate the disorder effects into the analysis. Since EXAFS is a locally sensitive method, a model containing several hundred or a few thousand atoms is usually sufficiently large to obtain accurate configuration-averaged EXAFS. If one is interested in just modelling the motion of atoms around their equilibrium positions, simulation of the behaviour of the system for a time period up to 100 ps is often sufficient (Okamoto, 2004; Kuzmin & Evarestov, 2009).

Figure 2 Flowchart of the classical molecular-dynamics method for EXAFS data analysis (Kuzmin & Evarestov, 2009).

As in the case of RMC, one also needs to provide an initial structure model for MD simulations. In the MD approach one also needs to set the velocities of the atoms. This is ensured by a so-called equilibration process, in which the velocities are assigned iteratively to each of the atoms in the model, so that, for instance, a distribution of velocities corresponds to the required temperature. Once thermodynamical equilibrium has been reached, the main run (production run) of the MD simulation can be started and the current atomic positions (snapshots) are recorded after given time intervals.

For each of the snapshots the corresponding theoretical EXAFS spectrum is calculated. Afterwards, all obtained configuration-averaged EXAFS spectra are averaged over time to obtain a single EXAFS spectrum, which includes contributions from thermal and static disorder, and therefore can be compared directly with the experimental EXAFS data (Palmer et al., 1996; Ferlat et al., 2005). An alternative approach for EXAFS calculations from MD data is to first calculate the partial RDFs for different atom pairs from the coordinates of MD simulations and then to integrate them to calculate the corresponding EXAFS spectrum (D'Angelo et al., 1994; Ferlat et al., 2005). While this approach is computationally less demanding, its applicability is limited to systems where multiple-scattering effects can be neglected.

If the agreement between experiment and theory is good, one can claim that the force-field model and selected structure model used reproduce the positions and dynamics of the atoms in the investigated material. Therefore, EXAFS can be used for the validation of force-field models (Kuzmin et al., 2016). An important point, especially for classical MD–EXAFS analysis, is that the information obtained on the disagreement between experimental and calculated EXAFS can also be used to modify the initial structure model or to optimize the force field (Anspoks et al., 2012). Thus, in principle, such classical MD analysis of EXAFS spectra can also be considered as an EXAFS fitting procedure, where the optimized variables are the parameters of the empirical force field and the structure parameters.

Some popular software codes for classical MD simulations include GULP (Gale & Rohl, 2003), LAMMPS (Plimpton, 1995) and DL_POLY (Todorov et al., 2006). For ab initio MD simulations, the CP2K (Lippert et al., 1999) and VASP (Kresse, 2000) codes, for example, can be employed.

Most RMC–EXAFS and MD–EXAFS studies are devoted to disordered materials (liquids, solutions, glasses and amorphous systems). Since in this case the information content of the EXAFS spectrum is often limited to the contribution of the first coordination shell only, EXAFS data alone do not allow one to obtain an unambiguous structure model, especially taking into account the broad and non-Gaussian distributions of bond lengths expected for such materials (Filipponi, 2001). Therefore, the information on atomic dynamics from classical MD simulations can significantly reduce the uncertainties of the analysis and this approach has successfully been applied, for example, to study ion solutions in water (D'Angelo et al., 1994; Palmer et al., 1996; Ferlat et al., 2005) and molten salts (Okamoto, 2004). Ab initio MD methods have also been applied for this purpose; for example, Car–Parrinello MD was applied for the interpretation of EXAFS data from Co²⁺ ions in water (Spezia et al., 2006), while in Ferlat et al. (2004) Born–Oppenheimer MD was found to be useful to determine the distributions of interatomic distances in In_100−xSe_x liquid alloys.

On the other hand, the ability of the RMC approach to combine the results of different experimental techniques is also of great importance in cases where the information available from EXAFS only is scarce. In RMC studies of disordered materials, EXAFS data are commonly used together with information from X-ray and neutron scattering experiments. Also, EXAFS data from several absorption edges are often used simultaneously. Examples of such combined RMC studies for disordered materials include work on chalcogenide materials (Pethes et al., 2016), metallic glasses (Kaban et al., 2014; Yang et al., 2014) and oxide glasses (Wicks et al., 1995).

For disordered materials the MD and RMC methods can be considered to be highly complementary. For instance, the results of ab initio MD–EXAFS simulations were compared with the results of RMC analysis in work on Ni–P metallic glasses (Sheng et al., 2006). Attempts to combine the MD and RMC methods have also been made. RMC can be used to generate an initial disordered structure model for ab initio MD simulations (Voleská et al., 2012). The existing knowledge of interatomic interactions, in turn, can be used as additional constraints in an RMC run to avoid nonphysical solutions (Gereben & Pusztai, 2012). This combined MD–RMC approach is also referred to as hybrid RMC (Opletal et al., 2002). Finally, RMC-type methods can be employed to refine an empirical interatomic potential (Soper, 1996; Bowron, 2008), which can be further used to generate a set of structural configurations by the Monte Carlo or MD procedure for the calculation of ensemble-averaged EXAFS, similarly to as is performed in the conventional MD–EXAFS approach.

Much richer information can be obtained from EXAFS spectra of crystalline materials, since in this case the distant coordination shells also contribute to the total EXAFS. From a technical point of view, simulation-based analysis of EXAFS data for crystalline materials is also more demanding, since the contributions of many strong multiple-scattering effects need to be included in the analysis. Nevertheless, the feasibility of MD–EXAFS and RMC–EXAFS analysis for contributions to the total EXAFS from distant atoms located at distances up to 10 Å from the absorbing atom has been demonstrated (Kuzmin et al., 2016; Timoshenko, Anspoks et al., 2014; Jonane et al., 2018).

As mentioned previously, the accuracy of even the most advanced EXAFS analysis for the determination of the equilibrium structure of a crystalline material will in most cases be lower than that provided by diffraction methods. Nevertheless, MD–EXAFS and RMC–EXAFS analysis can be a very powerful tool to investigate disorder phenomena in crystalline materials (Winterer et al., 2002; Németh et al., 2012; Jonane, Lazdins et al., 2016; Jonane, Timoshenko et al., 2016; Kalinko et al., 2016) and can shine some light, for example, on intriguing cases in which the displacements of atoms from their average, periodic structure result in unusual and technologically important properties of the material. As an example of such studies, one can mention the work of Krayzman, Levin and Tucker on perovskite-type materials (Levin et al., 2013), where the deviations of the local structure of perovskite from the crystallographic structure were investigated and correlated with the ferroelectric properties of the material using the RMC–EXAFS approach. Another example in which anharmonic effects and the anisotropy of atomic displacements play an important role is the so-called negative thermal expansion effect (NTE): contraction of the crystalline lattice upon temperature increase. A model material for such studies, rhenium trioxide, was investigated by means of the RMC–EXAFS method by Timoshenko, Kuzmin et al. (2014), and it was demonstrated that the RMC–EXAFS approach is indeed able to provide unique information on the correlations in atomic motion. The ab initio MD–EXAFS approach, in turn, appeared to be useful for studies of NTE in scandium fluoride (Bocharov et al., 2016).

Finally, the idea of using MD and RMC methods for the analysis of nanostructured materials seems to be very promising (Timoshenko et al., 2019). From the EXAFS point of view, many nano­materials can be considered as a crossover between completely disordered and crystalline structures. Broad and non-Gaussian distributions of bond lengths are expected in this case. The structure of such materials can also be quite different from the respective bulk structure due to structure-relaxation effects caused by the confined size and the presence of defects. The lack of long-range order limits the applicability of conventional structural analysis in this case (Billinge & Levin, 2007). At the same time, contributions from distant coordination shells to the total EXAFS and multiple-scattering effects can often still be important and may contain useful structural information. Hence, it is very appealing to apply simulation-based approaches to the analysis of EXAFS data in this case.

The first attempts to use RMC and EXAFS to refine the 3D structures of nanoparticles were the study of metallic Ni–Pt nanoparticles (Tupy et al., 2012), investigations of silver doping of CdSe nanocrystals (Kompch et al., 2015) and studies of cobalt and copper tungstate nanoclusters (Timoshenko et al., 2015) and of gold nanoparticles (Timoshenko & Frenkel, 2017).

Both the classical and ab initio MD approaches also appear to be useful for the EXAFS analysis of nano­structured materials. For example, a complex multi-step MD–EXAFS procedure was applied in Anspoks et al. (2012) to investigate the role of oxygen vacancies in the structure relaxation of nanosized NiO. Studies of metallic nanoparticles are another interesting example of the application of the MD–EXAFS method. While some attempts to use classical force-field models for the interpretation of EXAFS data from such systems are known (Clausen & Nørskov, 2000; Roscioni et al., 2011; Timoshenko & Frenkel, 2017; Timoshenko et al., 2017) and a good agreement between the results of experiment and simulation for the bare, unsupported nanoparticles can be obtained, the ab initio MD approach seems to be most promising in this case. The latter is able to account for the interactions between nanoparticles and surrounding ligands and various adsorbates, which may have a significant influence on the structure and functionality of nanoparticles. For instance, ab initio MD–EXAFS appeared to be useful for investigations of disorder within gold nanoparticles of 1–2 nm in size due to interactions with surrounding thiol ligands (Yancey et al., 2013; Chill et al., 2015).

Theoretical simulations can also be used for the analysis of XANES data. While the accuracy of ab initio XANES calculations is still limited and such calculations are also more time-consuming, examples of simulation-based interpretation of XANES spectra are known (Vila et al., 2008). For instance, the application of the classical MD method to the interpretation of XANES data from borosilicate glasses is discussed in Cabaret et al. (2001). At the same time, while the RMC-based analysis of XANES data is still prohibitive, XANES simulations for the final configuration obtained in RMC–EXAFS modelling can be used for additional validation of the obtained structural model (Sheng et al., 2006; Timoshenko et al., 2015; Jonane et al., 2019).

Numerical simulation of the 3D structure of materials for the interpretation of EXAFS spectra is a very powerful tool that allows one to exploit more of the information available from EXAFS data and to obtain a more detailed view of the structure of the material. It allows one to combine insights from different experimental techniques, and also to bridge the gap between the results of theoretical calculations and the experimental data. At the same time, there is still much room for further improvements of RMC-like and MD-like approaches. The development of fields such as evolutionary computation and machine learning (Timoshenko et al., 2018, 2020, 2023; Timoshenko & Frenkel, 2019) may contribute greatly to this process.