Chapter 6.25. XSpectra: a density-functional-theory-based plane-wave pseudopotential code for XANES calculation

XSpectra is a code dedicated to the calculation of X-ray absorption near-edge and pre-edge structures, i.e. the XANES region of an X-ray absorption spectrum. It is based on density-functional theory (DFT) and the Born–Oppenheimer approximation; it uses plane-wave basis sets, pseudopotentials and periodic boundary conditions.

XSpectra is distributed within the Quantum ESPRESSO suite of open-source codes (where ESPRESSO stands for opEn Source Package for Research in Electronic Structure) under the terms of the GNU General Public Licence (Giannozzi et al., 2009, 2017). It is used as a post-process for a self-consistent field (SCF) calculation.

Although DFT is a ground-state theory, it is used here to model excited states according to the final state rule, with the core-hole effects being taken into account within a supercell. More generally, the use of a supercell extends the applications to noncrystalline materials such as amorphous materials and finite systems. The standard implementation of DFT relies on an SCF calculation of the Kohn–Sham orbitals and energies using the exchange and correlation potential given in the local density approximation (LDA) or in the generalized gradient approximation (GGA) and including spin polarization if needed. This single-particle framework makes XSpectra well suited to the calculation of edges with delocalized final states, such as most of the K and L₁ edges and the L_2,3 and M_2,3 edges of elements belonging to the 4d and 5d series.

The first version of XSpectra was developed for K and L₁ edges by Gougoussis, Calandra, Seitsonen, Brouder et al. (2009), who generalized the norm-conserving pseudopotential method described by Taillefumier et al. (2002) to the case of ultrasoft pseudopotentials. The scheme was then extended to the calculation of L_2,3 edges by Bunău & Calandra (2013).

The fact that XSpectra belongs to the Quantum ESPRESSO integrated suite of codes has multiple advantages (Giannozzi et al., 2009). One of them is that a XANES simulation can be coupled with the calculation of other properties of materials under study within the same framework. The two examples given in Section 3 have been selected to illustrate (i) the combination of XANES with structure optimization and density-of-states (DOS) calculations and (ii) the interest in including quantum thermal fluctuations of nuclei.

In a single-particle DFT-based approach, the X-ray absorption cross section is given by where |ψ_i〉 is the single-electron initial state with energy E_i, |ψ_f〉 is the single-electron final (empty) state with energy E_f, is the incident X-ray beam energy and α is the fine-structure constant. In the electric quadrupole approximation, the transition operator is a sum of two terms, i.e. , where and are the electric dipole (E1) and electric quadrupole (E2) contributions, respectively. The quantities and k are the polarization-vector direction and the wavevector of the photon beam, respectively, and r is the position coordinate of the electron.

While the |ψ_i〉 initial states are simply described by core-level wavefunctions taken from an atomic ground-state calculation, the |ψ_f〉 empty states result from an SCF calculation for a supercell, which includes a core hole on the absorbing atom according to the final state rule. In a pseudopotential approach, the core electrons are considered with the nuclei as frozen ion cores, and only the chemically active valence electrons are explicitly treated in the calculations. Thus, the SCF run in Quantum ESPRESSO does not produce all-electron wavefunctions ψ (as a full-potential code does), but pseudo-wavefunctions . Since equation (1) requires all-electron final-state wavefunctions ψ_f, XSpectra performs an all-electron reconstruction, using the projector-augmented-wave (PAW) formalism described by Blöchl (1994). As shown by Taillefumier et al. (2002), the PAW formalism permits equation (1) to be rewritten as where are the eigenstates of the pseudo-Hamiltonian of the supercell. In equation (3), the wavefunctions are the all-electron partial waves centred on the R₀ absorbing atom position, the vectors form a complete set of projector functions and the n index refers to the angular momentum quantum numbers (l, m) and to an additional number ν, which is used if there is more than one projector per angular momentum channel. The PAW method requires pseudo-partial waves such that inside a spherical core region Ω_R, the so-called augmentation region. The vectors are then built in such a way that they are equal to zero outside Ω_R and satisfy the condition = δ_RR′δ_nn′. It is worth noting that equation (2) is obtained under the assumption that the initial state is localized on the absorbing atom, so that the overlap between and can be neglected if R ≠ R₀. Furthermore, even if equation (3) is in principle valid for a complete set of projectors, it converges after a few terms. In XSpectra, the all-electron partial waves are chosen as the solutions of the Schrödinger equation for the isolated atom and the PAW reconstruction is performed using its implementation in the GIPAW code of the Quantum ESPRESSO distribution, which was developed for nuclear magnetic resonance (NMR) and electronic paramagnetic resonance (EPR) calculations.

When dealing with supercells (typically 100 atoms), the calculation of XANES from equation (2) requires the calculation of many empty energy states ψ_f and thus many empty bands, a computationally expensive task that scales with , where N_G is the number of G vectors used in the simulation and N_f is the number of states up to the highest energy in the spectrum. To avoid this drawback, XSpectra builds a Lanczos basis using a powerful recurrence method (Haydock, 1980) that (i) transforms a Hermitian matrix into a tridiagonal form and (ii) permits the cross section to be rewritten as a continued fraction, such asIn equation (4), the operator S defined in the ultrasoft pseudopotential scheme reads where the q_R,nm integrated augmentation charges are given by q_R,nm = 〈φ_R,n|φ_R,m〉 − (Gougoussis, Calandra, Seitsonen & Mauri, 2009). The {aⁱ} and {bⁱ} sets of coefficients are the diagonal and subdiagonal terms of the matrix, respectively. It is worth noting that the use of norm-conserving pseudopotentials coincides with the S = 1 case. The number of terms of the continued fraction required for convergence strongly depends on the γ value, which is used as a broadening parameter that at least includes the core-hole lifetime. To improve the convergence, an analytical terminator is used to finish the continued fraction (Rocca et al., 2008).

With the Lanczos iterative technique, only the occupied bands need to be calculated, and the computational cost is reduced to , where N_L is the number of Lanczos iterations that are needed to converge. The computing time is considerably reduced compared with that required by the calculation of the final states ψ_f [equation (2)]. Furthermore, only two wavefunctions need to be stored.

The current version of XSpectra (5.4.0) calculates K, L₁ and L_2,3 edges using linear polarization of the X-ray beam (Gougoussis, Calandra, Seitsonen & Mauri, 2009; Bunău & Calandra, 2013). Both electric dipole (E1) and electric quadrupole (E2) transitions can be calculated (independently). Spin-polarized cross sections can be obtained for magnetic systems. Correlation effects can be simulated in a mean-field way using the Hubbard U correction (Cococcioni & de Gironcoli, 2005), which has been included in the code by Gougoussis, Calandra, Seitsonen, Brouder et al. (2009). Fig. 1 illustrates the type of theory–experiment agreement that may be expected using XSpectra in cases where a single-particle DFT approach is appropriate. XSpectra does not calculate absolute energies. By default, the zero of the calculated spectra is set to the Fermi level. Thus, comparison with experimental spectra needs a manual energy shift of the calculated spectrum. An energy-dependent broadening parameter, such as that described in section III of Bunău & Calandra (2013), can be used for γ in the continued fraction or in a convolution post-process of the cross section.

Figure 1 Experimental and calculated XANES spectra at the Si K edge in α-quartz (upper panel) and at the S L2,3 edges in the 2H polytype of MoS2 (lower panel). The hexagonal supercells used for the calculation are shown with the c axis vertical. The isotropic XANES spectrum, σiso, is given by (σ∥ + 2σ⊥)/3, where σ∥ and σ⊥ correspond to oriented along and normal to the c axis of the crystal, respectively. The X-ray natural linear dichroism (XNLD) signal is given by σ∥ − σ⊥. The quartz and MoS2 graphs were drawn from the data of Gougoussis, Calandra, Seitsonen & Mauri (2009) and Bunău & Calandra (2013), respectively.

Recent, ongoing and future developments of XSpectra are focused on (i) X-ray magnetic circular dichroism (XMCD) and X-ray natural circular dichroism (XNCD) at the K and L₁ edges (Bouldi et al., 2017), (ii) X-ray Raman spectroscopy (XRS; de Clermont Gallerande et al., 2018), (iii) electron-energy core-loss spectroscopy (ELNES) and (iv) the creation of an online GIPAW pseudopotential database dedicated to core-electron excitation spectroscopies.

Paramagnetic impurities are responsible for the colour of many allochromatic minerals. For instance, the substitution of Al³⁺ by Cr³⁺ induces the red colour of ruby (α-Al₂O₃) and spinel (MgAl₂O₄) and the green colour of emerald (Be₃Al₂Si₆O₁₈). The colouring mechanism is intimately connected to the local structure of Cr³⁺, which cannot be assessed by X-ray diffraction. Thus, use of a local probe such as X-ray absorption is mandatory.

The structural and electronic environment of Cr³⁺ in spinel has been investigated by combining EXAFS and XANES at the Cr K edge with first-principles calculations (Juhin et al., 2007, 2008). The atomic positions of a 2 × 2 × 2 rhombohedral supercell (1 Cr, 31 Al, 16 Mg, 64 O) were relaxed by total energy minimization. The full relaxation of the coordination sphere is observed with a theoretical Cr—O distance close to that in MgCr₂O₄ and in good agreement with the EXAFS analysis. Despite some angular site distortion, the point group of the host site is conserved. Beyond the first shell, the Cr for Al substitution slightly increases the distance of the Al second shell. The XANES spectrum is found to be sensitive to these structural relaxation effects, as shown in Fig. 2, which compares the calculated spectra performed with the relaxed supercell (solid line) and with the nonrelaxed supercell (dashed line) with the experimental spectra. The relaxation does significantly impact the A, B, C and E spectral features, leading to a better agreement with experiment. To conclude, this theoretical and experimental coupled study proves, at least to a certain extent, the reliability of the MgAl₂O₄:Cr³⁺ relaxed structural model.

Figure 2 Sensitivity of XANES on local relaxation around a Cr impurity in red spinel (MgAl2O4:Cr3+). In the nonrelaxed calculation, Cr substitutes for an Al in its exact position. Data were extracted from Juhin et al. (2007).

The Cr K pre-edge region (framed in Fig. 2) contains information about the Cr 3d states. Fig. 3 is focused on this pre-edge and shows the E2 origin of the two P₁ and P₂ features. Indeed, only pure 1s→3d transitions are expected in this case since (i) the absorbing atom site is centrosymmetric (no in-site E1 transition owing to p–d hybridization) and (ii) Cr is the unique 3d element isolated in the supercell (no off-site E1 transition owing to p–d mixing with neighbours). The spin-polarized 3d DOS projected for the absorbing Cr (lower panel) allows the interpretation of the pre-edge in terms of monoelectronic transitions. In the single-particle picture of Cr³⁺ in an octahedral site, the Cr 3d t_2g states are occupied by three electrons, while the e_g states are empty. This description coincides with the DOS calculations: P₁ originates from both 1s→3d and 1s→3d transitions and P₂ is explained by 1s→3d transitions. Besides, this example shows that a DFT single-particle framework enables the interpretation of K pre-edge features of transition-metal elements. More examples are given in Cabaret et al. (2010), where the limitations of a single-particle DFT-based framework to calculate transition-metal K pre-edges are widely discussed. In particular, E2 and local E1 transitions were systematically found at too high an energy with respect to the edge. This kind of drawback, which is related to the modelling of the 1s core hole–electron interaction in a supercell, could be avoided using Bethe–Salpeter-based approaches (Gilmore et al., 2015).

Figure 3 Cr K pre-edge analysis in MgAl2O4:Cr3+ (). Top panel: decomposition of the total calculated spectrum into the E1 and E2 components. Data were extracted from Juhin et al. (2008). Bottom panel: spin-polarized 3d density of states projected for the Cr impurity of the supercell. The same broadening parameter was used for the XANES and DOS calculations, erasing the electronic gap in the DOS panel.

Most XANES calculations in crystals consider the atoms to be fixed at their equilibrium positions as determined by X-ray diffraction. However, at least two arguments indicate that nuclear motion may have a significant effect on XANES, notably in low-Z cation oxides: (i) the usual XANES calculation fails to reproduce some pre-edge features at the cation K edge and (ii) temperature induces significant spectral modifications, especially in the pre-edge (Manuel et al., 2012). It has been shown that such pre-edge peaks are related to a combined effect of the core hole–electron interaction and 1s→3d transitions that become allowed in the electric dipole approximation through s–p hybridization induced by vibrations (Manuel et al., 2012; Nemausat et al., 2015).

Nemausat et al. (2015) developed a theoretical framework to include the thermal quantum fluctuations of nuclei in XANES calculations. From phonon calculations performed in the quasi-harmonic approximation (QHA) using the PHonon code of Quantum ESPRESSO, non-equilibrium atomic configurations are generated such that they obey quantum statistics at finite temperature (Errea et al., 2013). For each configuration, an individual spectrum is calculated. Then the temperature-dependent theoretical XANES spectrum is obtained by averaging the individual spectra that have been core-level shifted beforehand. The case of the Al K edge in corundum (α-Al₂O₃) at 300 K is presented in Fig. 4 (σ_∥ orientation), where all the 90 configuration spectra ensuring convergence are displayed in the background. The QHA spectrum (solid red line) globally is in better agreement with experiment than the standard equilibrium spectrum (black solid line). This is especially true in the pre-edge region, as shown in the inset, which emphasizes the absence of any P feature if vibrations are not taken into account. Similar calculations performed at 0 K also exhibit a pre-edge peak (Nemausat et al., 2017). This highlights the quantum nature of the nuclear motion at the origin of peak P.

Figure 4 The impact of quasi-harmonic phonons on the Al K-edge XANES spectrum of α-Al2O3. The σ∥ experimental spectrum is compared with two calculated spectra obtained either with the atoms fixed at their equilibrium positions (black) or in the QHA framework (red) according to Nemausat et al. (2015) from an average of the configuration spectra shown in the background (grey lines). Inset: enlargement of the pre-edge region, highlighting the vibrational origin of the pre-edge peak P.