International
Tables for Crystallography Volume I X-ray absorption spectroscopy and related techniques Edited by C. T. Chantler, F. Boscherini and B. Bunker © International Union of Crystallography 2020 |
International Tables for Crystallography (2020). Vol. I. Early view chapter
https://doi.org/10.1107/S1574870720003328 XSpectra : a density-functional-theory-based plane-wave pseudopotential code for XANES calculation
a
Sorbonne Université, Campus Pierre et Marie Curie, Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie, UMR7590, Case Courrier 115, 4 Place Jussieu, 75005 Paris, France, and bSorbonne Université, CNRS, Institut des Nanosciences de Paris, UMR7588, 75005 Paris, France The XSpectra code calculates XANES spectra at the K, L1 and L2,3 edges, including both electric dipole and quadrupole transitions. Density-functional theory and the Born–Oppenheimer approximation are used. The implementation relies on plane-wave basis sets, pseudopotentials and periodic boundary conditions. XSpectra belongs to the Quantum ESPRESSO distribution, an integrated suite of open-source computer codes for electronic structure calculations. XSpectra reads the self-consistent charge density produced by the PWscf code of the Quantum ESPRESSO distribution and acts as a post-processing tool. Core-hole effects are included by considering a supercell containing one absorbing atom, the pseudopotential of which is generated with a core hole in the appropriate core level. The all-electron final state wavefunctions are constructed using the projector augmented-wave method. The XANES cross section is calculated using the Lanczos method and a continued fraction expansion, avoiding the explicit calculation of empty states. The capabilities of the code are illustrated via selected examples. Keywords: XANES; XSpectra; first-principles calculations; density-functional theory; DFT; plane waves; pseudopotentials; Lanczos method. |
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 L1 edges and the L2,3 and M2,3 edges of elements belonging to the 4d and 5d series.
The first version of XSpectra was developed for K and L1 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 L2,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 Ei, |ψf〉 is the single-electron final (empty) state with energy Ef, 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 R0 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 ≠ R0. 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 NG is the number of G vectors used in the simulation and Nf 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 qR,nm integrated augmentation charges are given by qR,nm = 〈φR,n|φR,m〉 − (Gougoussis, Calandra, Seitsonen & Mauri, 2009). The {ai} and {bi} 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 NL 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, L1 and L2,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.
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 L1 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.
The theoretical XANES spectrum calculation is a two-step process. Firstly, an SCF run is performed using PWscf for a supercell including the core hole on the absorbing atom. Secondly, the thus-generated SCF charge density is used to calculate the XANES cross section with XSpectra. An XSpectra run is performed considering either E1 transitions for a given orientation of in the supercell or E2 transitions requiring both and k orientations provided .
The core-hole electron interaction is usually modelled using the full core-hole approach (FCH): the one-electron charge deficit due to the core hole is compensated by a background charge uniformly distributed within the supercell volume. Alternatively, half-core-hole (HCH) and excited-core-hole (XCH) approaches are possible. The latter consists of adding one electron localized on the first empty band, leading to a core-hole screening stronger than in the FCH approach. XCH was compared with FCH and HCH at the Li, B and C K edges of XRS spectra in the case of Li-bearing reference compounds and was found to slightly improve the agreement with experiment (de Clermont Gallerande et al., 2018). All of the theoretical spectra presented here are FCH calculations.
The main calculation parameters to be converged are the following: (i) the cutoff energy for the plane-wave basis set, (ii) the k-point grid to sample the first Brillouin zone (BZ) for the SCF run, (iii) the k-point grid for the XSpectra calculation and (iv) the supercell size. Parameters (i), (ii) and (iii) can be converged on the unit cell. The supercell size has to be large enough to provide a sufficient spatial separation (typically 8 Å) between the absorbing atom and its periodic image. The k-point mesh to be used in the supercell in (ii) and (iii) is then obtained by keeping the density of k-points in the supercell BZ identical to that in the unit cell.
For each k-point of the BZ, XSpectra calculates [equation (3)] and uses it to perform the Lanczos algorithm. The convergence of the Lanczos basis set, which strongly depends on the γ value, is ensured at each k-point. The k-point loop is performed for each initial state, i.e. one for the K and L1 edges, two for the L2 edge and four for the L3 edge. Then, the XANES spectrum to be compared with experiment is obtained using equation (4) by gathering the contributions of all k-points.
In the case of more than one non-equivalent atomic site for the absorbing atom, the spectrum is the weighted average of individual calculated spectra according to the mutiplicity of each absorbing atom site. Since the energy zeros of individual spectra are not identical, a core-level shift has to be evaluated before performing the individual spectra average. This can be perfomed by using the method described in Schwartz et al. (2009).
The pseudopotential files should include the φR,n partial waves and pseudo-partial waves needed to calculate the PAW projectors. In general, two linearly independent projectors by the (l, m) channel are sufficient to ensure the spectrum convergence. Bunău & Calandra (2013) show that three projectors are by far sufficient to describe the full L2 spectrum of Cu 50 eV above the edge. Norm-conserving and ultrasoft pseudopotentials produce the same results but with substantially reduced cutoff energy using the ultrasoft pseudopotentials. No clear evidence of the impact on the absorption cross section of the use of either LDA or GGA exchange-correlation functionals has yet been established.
Paramagnetic impurities are responsible for the colour of many allochromatic minerals. For instance, the substitution of Al3+ by Cr3+ induces the red colour of ruby (α-Al2O3) and spinel (MgAl2O4) and the green colour of emerald (Be3Al2Si6O18). The colouring mechanism is intimately connected to the local structure of Cr3+, 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 Cr3+ 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 MgCr2O4 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 MgAl2O4:Cr3+ relaxed structural model.
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 P1 and P2 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 Cr3+ in an octahedral site, the Cr 3d t2g states are occupied by three electrons, while the eg states are empty. This description coincides with the DOS calculations: P1 originates from both 1s→3d and 1s→3d transitions and P2 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).
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 (α-Al2O3) 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.
Acknowledgements
We greatly acknowledge the main contributors to the XSpectra code: Christos Gougoussis, Ari P. Seitsonen and Oana Bunau. Special thanks go to Nadejda Bouldi, Nora Vollmers, Christian Brouder and Uwe Gerstmann, the developers of the XMCD part, and to Guillaume Radtke, Emmanuelle de Clermont Gallerande, Steven Delhommaye and Lorenzo Paulatto, the developers of the XRS part.
References
Blöchl, P. (1994). Phys. Rev. B, 50, 17953–17979.Google ScholarBouldi, N., Vollmers, N. J., Delpy-Laplanche, C. G., Joly, Y., Juhin, A., Sainctavit, P., Brouder, C., Calandra, M., Paulatto, L., Mauri, F. & Gerstmann, U. (2017). Phys. Rev. B, 96, 085123.Google Scholar
Bunău, O. & Calandra, M. (2013). Phys. Rev. B, 87, 205105.Google Scholar
Cabaret, D., Bordage, A., Juhin, A., Arfaoui, M. & Gaudry, E. (2010). Phys. Chem. Chem. Phys. 12, 5619–5633.Google Scholar
Clermont Gallerande, E. de, Cabaret, D., Lelong, G., Brouder, C., Attaiaa, M.-B., Paulatto, L., Gilmore, K., Sahle, C. J. & Radtke, G. (2018). Phys. Rev. B, 98, 214104.Google Scholar
Cococcioni, M. & de Gironcoli, S. (2005). Phys. Rev. B, 71, 035105.Google Scholar
Errea, I., Calandra, M. & Mauri, F. (2013). Phys. Rev. Lett. 111, 177002.Google Scholar
Giannozzi, P., Andreussi, A., Brumme, T., Bunau, O., Buongiorno Nardelli, M., Calandra, M., Car, R., Cavazzoni, C., Ceresoli, D., Cococcioni, M., Colonna, N., Carnimeo, I., Dal Corso, A., de Gironcoli, S., Delugas, P., DiStasio, R. A., Ferretti, A., Floris, A., Fratesi, G., Fugallo, G., Gebauer, R., Gerstmann, U., Giustino, F., Gorni, T., Jia, J., Kawamura, M., Ko, H.-Y., Kokalj, A., Küçükbenli, E., Lazzeri, M., Marsili, M., Marzari, N., Mauri, F., Nguyen, N. L., Nguyen, H.-V., Otero-de-la-Roza, A., Paulatto, L., Poncé, S., Rocca, D., Sabatini, R., Santra, B., Schlipf, M., Seitsonen, A. P., Smogunov, A., Timrov, I., Thonhauser, T., Umari, P., Vast, N., Wu, X. & Baroni, S. (2017). J. Phys. Condens. Matter, 29, 465901.Google Scholar
Giannozzi, P., Baroni, S., Bonini, N., Calandra, M., Car, R., Cavazzoni, C., Ceresoli, D., Chiarotti, G. L., Cococcioni, M., Dabo, I., Dal Corso, A., de Gironcoli, S., Fabris, S., Fratesi, G., Gebauer, R., Gerstmann, U., Gougoussis, C., Kokalj, A., Lazzeri, M., Martin-Samos, L., Marzari, N., Mauri, F., Mazzarello, R., Paolini, S., Pasquarello, A., Paulatto, L., Sbraccia, C., Scandolo, S., Sclauzero, G., Seitsonen, A. P., Smogunov, A., Umari, P. & Wentzcovitch, R. M. (2009). J. Phys. Condens. Matter, 21, 395502.Google Scholar
Gilmore, K., Vinson, J., Shirley, E., Prendergast, D., Pemmaraju, C., Kas, J., Vila, F. & Rehr, J. (2015). Comput. Phys. Commun. 197, 109–117.Google Scholar
Gougoussis, C., Calandra, M., Seitsonen, A. P., Brouder, C., Shukla, A. & Mauri, F. (2009). Phys. Rev. B, 79, 045118.Google Scholar
Gougoussis, C., Calandra, M., Seitsonen, A. P. & Mauri, F. (2009). Phys. Rev. B, 80, 075102.Google Scholar
Haydock, R. (1980). Comput. Phys. Commun. 20, 11–16.Google Scholar
Juhin, A., Brouder, C., Arrio, M.-A., Cabaret, D., Sainctavit, P., Balan, E., Bordage, A., Seitsonen, A. P., Calas, G., Eeckhout, S. G. & Glatzel, P. (2008). Phys. Rev. B, 78, 195103.Google Scholar
Juhin, A., Calas, G., Cabaret, D., Galoisy, L. & Hazemann, J.-L. (2007). Phys. Rev. B, 76, 054105.Google Scholar
Manuel, D., Cabaret, D., Brouder, C., Sainctavit, P., Bordage, A. & Trcera, N. (2012). Phys. Rev. B, 85, 224108.Google Scholar
Nemausat, R., Cabaret, D., Gervais, C., Brouder, C., Trcera, N., Bordage, A., Errea, I. & Mauri, F. (2015). Phys. Rev. B, 92, 144310.Google Scholar
Nemausat, R., Gervais, C., Brouder, C., Trcera, N., Bordage, A., Coelho-Diogo, C., Florian, P., Rakhmatullin, A., Errea, I., Paulatto, L., Lazzeri, M. & Cabaret, D. (2017). Phys. Chem. Chem. Phys. 19, 6246–6256.Google Scholar
Rocca, D., Gebauer, R. Y., Saad, Y. & Baroni, S. (2008). J. Chem. Phys. 128, 154105.Google Scholar
Schwartz, C. P., Uejio, J. S., Saykally, R. J. & Prendergast, D. (2009). J. Chem. Phys. 130, 184109.Google Scholar
Taillefumier, M., Cabaret, D., Flank, A.-M. & Mauri, F. (2002). Phys. Rev. B, 66, 195107.Google Scholar