Reduced-Shifted Conjugate-Gradient Method with Seed Switching for Calculating X-ray Absorption Spectra

Abstract

We propose a new algorithm for numerical evaluation of x-ray absorption spectra based on the Green's function formalism in this paper. The method stands on the reduced-shifted conjugate-gradient (RSCG) method, which enables us to obtain the Green's functions with different energy-shifts simultaneously. The seed switching technique is combined into the RSCG method to improve the numerical accuracy and robustness of the algorithm, which we refer to the RSCG-SS method. The RSCG-SS method is applied for calculations of x-ray absorption spectra at the L_2,3-edges of transition metals in simple oxides as a benchmark. The theoretical spectra obtained by the new algorithms are identical to those obtained by the full-diagonalization method following the Fermi's golden rule. The RSCG-SS method can also be applied for the calculation of spectra with a sizable Hamiltonian matrix (∼300,000) with reasonable computational costs.

1. Introduction

The x-ray absorption spectroscopy (XAS) monitors electronic excitation of a core electron to unoccupied states via photon absorption. XAS becomes an indispensable tool in modern science for the characterization of the local electronic and atomic structures in materials.^1,2) Various theoretical approaches with different levels of approximation for the handling of electronic exchange-correlation interactions have been developed to reproduce and interpret measured XAS spectra from electronic structure calculations. One-particle calculations based on the multiple scattering theory or density functional theory (DFT) calculations, including core-hole effects, work well when the interactions between electrons and core-hole are weak.^3–12) Most of the K-edge and L_2,3-edges XAS of typical elements have been quantitatively reproduced by these methods. At the low-energy absorption edge below 100 eV, the interaction between the core-hole and the excited electron, namely the exciton, affect the XAS spectral shape. The two-particle approaches involving the solution of the Bethe-Salpeter equation (BSE) that can explicitly handle the excitonic effects, should be utilized.^13–16) In the case of L_2,3-edges of transition metals (TMs) or M_4,5-edges of rare-earth (RE) elements, which are dominated by the excitation into spatially localized 3d and 4f levels, the spectral shapes are significantly modified from the TM-3d and RE-4f density of states; this is known as the multiplet effects. The many-particle method that can explicitly deal with the exchange-correlation interaction among core-hole and 3d and 4f electrons are necessary.^17,18) The configuration interaction (CI) method and dynamically mean-field theory (DMFT) have successfully reproduced those spectra.^19–24) We refer to the reviews for a more detailed comparison of these theoretical approaches.^25–27)

In this work, the ab-initio CI method for simulating the multiplet structures at the L_2,3-edges of TM compounds,¹⁹⁾ which was developed by the present authors, was reconsidered. Theoretically, the absorption coefficient μ(ω) at the incident x-ray energy ħω is proportional to the photo-absorption cross-section (PACS) σ_abs(ω). According to the Fermi’s golden rule, σ_abs(ω) can be written in SI unit as,   

\begin{equation} \sigma_{\text{abs}}(\omega) = \sum_{f}4\pi^{2}\alpha |\langle \Psi_{f} | \skew3\hat{A} | \Psi_{i} \rangle |^{2}\delta (E_{f} - E_{i} - \hbar \omega), \end{equation} (1)

where, α is the fine-structure constant. Ψ_i and Ψ_f are the many-electron wavefunctions corresponding to the initial and final states, respectively, while E_i and E_f are their energies, respectively. $\skew3\hat{A}$ is the transition operator, which represents the interaction of photons with the electrons. In the original ab-initio CI method, many-electron wavefunctions Ψ_i and Ψ_f are expanded by Slater determinants constituting the relativistic molecular orbitals with different electronic configurations. The many-electron Hamiltonian matrix is constructed in terms of Slater determinants and fully diagonalized to obtain the initial and final states of the XAS excitation process. Then, PACS is evaluated directly by applying eq. (1). However, this approach is feasible only when the dimension of Hamiltonian matrix N is kept small as the computational time and required memory space scale as O(N³) and O(N²), respectively.

In the previous paper, we reported the CI method using the iterative algorithms to handle the larger matrix size.²⁸⁾ In this approach, PACS is evaluated solely from the initial states using the spectral function, which is the imaginary part of the Green's function:   

\begin{equation} \sigma_{\text{abs}}(\omega) = \lim_{\eta \to 0}\frac{4\pi^{2}\alpha }{\omega } \text{Im}\,G_{A}(\omega ,\eta), \end{equation} (2)

  

\begin{equation} G_{A}(\omega ,\eta) = - \frac{1}{\pi }\langle \Psi_{i} | \skew3\hat{A}^{\dagger }(zI - H)^{ - 1}\skew3\hat{A} | \Psi_{i}\rangle, \end{equation} (3)

where H is the Hamiltonian matrix, I is the identity matrix of dimension N, and z = ħω + E_i + iη. During the evaluation of eq. (2), linear equations (zI − H) x = b are solved for a set of different diagonal shift zI where the right-hand-side vector $\boldsymbol{{b}} = \skew3\hat{A}| \Psi _{i} \rangle $ is common to all equations. Then, the value of the spectral function is calculated as a dot product of two vectors G_A = b^†x. We adopted the generalized Davidson method for obtaining initial states, (E_i, Ψ_i), and the Lanczos method for evaluating the G_A(ω, η). Similar approaches were adopted for XAS simulation within the framework of multiple scattering theory,²⁹⁾ and plane-wave basis DFT calculations.³⁰⁾ Though the Lanczos method works well for many cases, care must be taken that it is prone to numerical instability in the finite precision arithmetic. A more robust algorithm would be preferable.

Recently, Nagai et al., proposed the reduced-shifted conjugate-gradient (RSCG) method, which can compute Green's function of the form G_αβ = ⟨α|(zI − H)⁻¹|β⟩ for multiple shift values simultaneously.³¹⁾ In the RSCG method, values of Green's function are computed directly without explicitly obtaining solution vectors of linear equations (zI − H) x = b. Hence, the memory consumption is significantly reduced. Yamamoto et al., introduced “the seed switching” into the shifted conjugate-orthogonal-conjugate-gradient (COCG) method to accelerate the convergence of the solutions of shifted linear equations for a real Hamiltonian matrix with complex shifts.³²⁾ There, the seed switching plays a vital role to improve the numerical robustness and to ensure the convergence of Green's function at every energy region of interest.

In the present work, we developed a new optimized algorithm for evaluating a spectrum function in eq. (3) by combining the RSCG method with seed switching (RSCG-SS). We applied this RSCG-SS method for the calculations of TM-L_2,3 XAS spectra of simple TM oxides as benchmarks. Theoretical spectra obtained by the proposed algorithm are found to be identical to those obtained by the full-diagonalization method. The proposed method can compute XAS spectra in the desired accuracy, even if the Hamiltonian matrix becomes large (∼300,000).

2. Theory

Let us consider a series of linear equations for matrices that differ only in diagonal components, so-called shifted linear systems,   

\begin{equation} (\sigma I + A)\boldsymbol{{x}}(\sigma) = \boldsymbol{{b}} \end{equation} (4)

where $A \in \mathbb{C}^{N \times N}$, $\boldsymbol{{x}}( \sigma ) \in \mathbb{C}^{N}$, $\boldsymbol{{b}} \in \mathbb{C}^{N}$, and $\sigma \in \mathbb{C}$. We assume that A is regular and a Hermitian matrix. We look for approximate solutions of eq. (4) in the Krylov subspace: $\mathcal{K}_{k}( A,\boldsymbol{{v}} )\ {:}{=}\ \text{span}\{ \boldsymbol{{v}},A\boldsymbol{{v}}, \ldots ,A^{k - 1}\boldsymbol{{v}} \}$.³³⁾ By using the fact that Krylov subspace is invariant under the scalar shift of diagonal elements with fixed vector, i.e., $\mathcal{K}_{k}( \sigma I + A,\boldsymbol{{v}} ) = \mathcal{K}_{k}( A,\boldsymbol{{v}} )$, we can reuse the Krylov subspace created at the fixed scalar shift, σ = σ^ref, to compute the approximate solutions for other σ. In ref. 31), the shifted CG method is applied to solve these shifted linear systems. Since the CG method is only applicable when σI + A is Hermitian, σ_ref should be taken on the real axis. The CG algorithm is constructed so that the residual vector r_k = b − Ax_k is perpendicular to the Krylov subspace: $\boldsymbol{{r}}_{k} \,\bot\, \mathcal{K}_{k}( A,\boldsymbol{{v}} )$. As the Krylov subspace is invariant for arbitrary shifts, the residual vectors of shifted linear systems are colinear to that of the seed system: $\boldsymbol{{r}}_{k}( \sigma ) = \rho _{k}( \sigma )\boldsymbol{{r}}_{k}^{\text{ref}}$.

For evaluation of Green’s function in eq. (3), it is not necessary to hold x(σ) explicitly: what we want to calculate is the values of g(σ_j) = b^†x(σ_j) with distinct scalar shifts $\{ \sigma _{j};j = 1, \ldots ,M\} $. The RSCG algorithm³¹⁾ is obtained by reformulating the equations in the shifted CG algorithm in terms of g(σ_j) rather than x(σ_j).

Here, we describe the RSCG-SS method in detail. As the first step of the RSCG, the linear equation in eq. (4) is solved at the reference scalar shift, $\sigma = \sigma _{\text{ref}} \in \mathbb{R}$. Hereafter, we refer to this equation as “the seed system.” The seed system is iteratively solved as follows under the initial conditions $\alpha _{0}^{\text{ref}} = 1$, $\beta _{0}^{\text{ref}} = 0$, $\boldsymbol{{r}}_{0}^{\text{ref}} = \boldsymbol{{b}}$:   

\begin{equation} \alpha_{k}^{\text{ref}} = \left(\frac{\beta_{k - 1}^{\text{ref}}}{\alpha_{k - 1}^{\text{ref}}} + \frac{(\boldsymbol{{r}}_{k}^{\text{ref}})^{\dagger }(\sigma_{\text{ref}}I + A)\boldsymbol{{r}}_{k}^{\text{ref}}}{(\boldsymbol{{r}}_{k}^{\text{ref}})^{\dagger }\boldsymbol{{r}}_{k}^{\text{ref}}}\right)^{ - 1}, \end{equation} (5)

  

\begin{equation} \theta_{k}^{\text{ref}} = \alpha_{k}^{\text{ref}}\frac{\beta_{k - 1}^{\text{ref}}}{\alpha_{k - 1}^{\text{ref}}}, \end{equation} (6)

  

\begin{equation} \boldsymbol{{r}}_{k + 1}^{\text{ref}} = (1 + \theta_{k}^{\text{ref}})\boldsymbol{{r}}_{k}^{\text{ref}} - \theta_{k}^{\text{ref}}\boldsymbol{{r}}_{k}^{\text{ref}} - \alpha_{k}^{\text{ref}}A\boldsymbol{{r}}_{k}^{\text{ref}}, \end{equation} (7)

  

\begin{equation} \beta_{k}^{\text{ref}} = \frac{(\boldsymbol{{r}}_{k + 1}^{\text{ref}})^{\dagger }\boldsymbol{{r}}_{k + 1}^{\text{ref}}}{(\boldsymbol{{r}}_{k}^{\text{ref}})^{\dagger }\boldsymbol{{r}}_{k}^{\text{ref}}}. \end{equation} (8)

Note that the solution vector $\boldsymbol{{x}}_{k + 1}^{\text{ref}}$ is not explicitly calculated: only the residual vector $\boldsymbol{{r}}_{k + 1}^{\text{ref}}$ is required for calculating Green’s function in the latter step. The equations for shifted linear systems are constructed in terms of the values of Green’s functions, g(σ₁), …, g(σ_M), are expressed as follows:   

\begin{align} &\rho_{k + 1}(\sigma_{j}) \\ &= \frac{\alpha_{k - 1}^{\text{ref}}\rho_{k}(\sigma_{j})\rho_{k - 1}(\sigma_{j})}{\alpha_{k - 1}^{\text{ref}}\rho_{k - 1}(\sigma_{j})(1 + \alpha_{k}^{\text{ref}}\sigma_{j}) + \alpha_{k}^{\text{ref}}\beta_{k - 1}^{\text{ref}}(\rho_{k - 1}(\sigma_{j}) - \rho_{k}(\sigma_{j}))}, \end{align} (9)

  

\begin{equation} \text{g}_{k + 1}(\sigma_{j}) = \alpha_{k}^{\text{ref}}\frac{\rho_{k + 1}(\sigma_{j})}{\rho_{k}(\sigma_{j})}\text{q}_{k}(\sigma_{j}), \end{equation} (10)

  

\begin{equation} \text{q}_{k + 1}(\sigma_{j}) = \beta_{k}^{\text{ref}}\left(\frac{\rho_{k + 1}(\sigma_{j})}{\rho_{k}(\sigma_{j})}\right)^{2}\text{q}_{k}(\sigma_{j}), \end{equation} (11)

with ρ₋₁(σ_j) = ρ₀(σ_j) = 1, g₀(σ_j) = 0, and q₀(σ_j) = b^†b for j = 1, …, M. The convergence of g_k+1(σ_j) is checked by monitoring the norms of residual vectors of shifted linear systems, which can be evaluated as,   

\begin{equation} r_{k + 1}(\sigma_{j}) = \| \boldsymbol{{r}}_{k + 1}(\sigma_{j})\|_{2} = | \rho_{k + 1}(\sigma_{j}) |\| \boldsymbol{{r}}_{k + 1}^{\text{ref}} \|_{2}. \end{equation} (12)

It should be noted that the matrix-vector product, which is the most time-consuming part of this algorithm for large matrix, is required only once per iteration for updating $\boldsymbol{{r}}_{k}^{\text{ref}}$ in the reference system. All the quantities in shifted linear systems, including the values of the Green’s function, can be obtained by simple scalar arithmetic.

The convergence behavior of the RSCG algorithm described above strongly depends on the choice of σ_ref. The convergence rate of spectrum function g_k+1(σ_j) is significantly reduced or is diverged in the worst case, when the residual norm of the seed system, $\| \boldsymbol{{r}}_{k + 1}^{\text{ref}} \|_{2}$, becomes small during the RSCG iterations. The seed switching technique is embedded to avoid such a situation, in which we replace the reference scalar shift, σ^ref, to a new one, $\sigma _{\text{new}}^{\text{ref}}$, when $\| \boldsymbol{{r}}_{k + 1}^{\text{ref}} \|_{2}$ becomes smaller than the prescribed threshold value, 0 < ϵ ≪ 1. To determine appropriate $\sigma _{\text{ref}}^{\text{new}}$ automatically, we introduce a set of auxiliary shifted systems, $( \sigma _{j}^{\text{aux}}I + A )\boldsymbol{{x}}^{\text{aux}}( \sigma _{j}^{\text{aux}} ) = \boldsymbol{{b}}$, with distinct real shifts, $\{ \sigma _{j}^{\text{aux}} \in \mathbb{R};j = 1, \ldots ,M' \} $, and track the scalar factor $\rho _{k}^{\text{aux}}( \sigma _{j}^{\text{aux}} )$ and residual norm, $r_{k}^{\text{aux}}( \sigma _{j}^{\text{aux}} ) = \| b - ( \sigma _{j}^{\text{aux}}I + A )\boldsymbol{{x}}_{k}^{\text{aux}}( \sigma _{j}^{\text{aux}} ) \|_{2}$, in each iteration:   

\begin{equation} \rho_{k + 1}^{\text{aux}}(\sigma_{j}^{\text{aux}}) = \frac{\alpha_{k - 1}^{\text{ref}}\rho_{k}^{\text{aux}}(\sigma_{j}^{\text{aux}})\rho_{k - 1}^{\text{aux}}(\sigma_{j}^{\text{aux}})}{\alpha_{k - 1}^{\text{ref}}\rho_{k - 1}^{\text{aux}}(\sigma_{j}^{\text{aux}})(1 + \alpha_{k}^{\text{ref}}\sigma_{j}^{\text{aux}}) + \alpha_{k}^{\text{ref}}\beta_{k - 1}^{\text{ref}}(\rho_{k - 1}^{\text{aux}}(\sigma_{j}^{\text{aux}}) - \rho_{k}^{\text{aux}}(\sigma_{j}^{\text{aux}}))}, \end{equation} (13)

  

\begin{equation} r_{k + 1}^{\text{aux}}(\sigma_{j}^{\text{aux}}) = | \rho_{k + 1}^{\text{aux}}(\sigma_{j}^{\text{aux}}) |\| \boldsymbol{{r}}_{k + 1}^{\text{ref}} \|_{2}. \end{equation} (14)

The new shift is set to be $\sigma _{\text{new}}^{\text{ref}} = \sigma _{s}^{\text{aux}}$ with $s = \arg \underset{j}{\max} r_{k + 1}^{\text{aux}}(\sigma _{j}^{\text{aux}} )$. During the seed switching step, internal parameters in the seed system and shifted linear systems are rescaled as follows:   

\begin{equation} \alpha_{k}^{\text{ref}} \leftarrow \frac{\rho_{k + 1}^{\text{aux}}(\sigma_{\text{new}}^{\text{ref}})}{\rho_{k}^{\text{aux}}(\sigma_{\text{new}}^{\text{ref}})}\alpha_{k}^{\text{ref}}, \end{equation} (15)

  

\begin{equation} \beta_{k}^{\text{ref}} \leftarrow \left(\frac{\rho_{k + 1}^{\text{aux}}(\sigma_{\text{new}}^{\text{ref}})}{\rho_{k}^{\text{aux}}(\sigma_{\text{new}}^{\text{ref}})}\right)^{2}\beta_{k}^{\text{ref}}, \end{equation} (16)

  

\begin{equation} \boldsymbol{{r}}_{k}^{\text{ref}} \leftarrow \rho_{k}^{\text{aux}}(\sigma_{\text{new}}^{\text{ref}})\boldsymbol{{r}}_{k}^{\text{ref}}, \end{equation} (17)

  

\begin{equation} \rho_{k}(\sigma_{j}) \leftarrow \frac{\rho_{k}(\sigma_{j})}{\rho_{k}^{\text{aux}}(\sigma_{\text{new}}^{\text{ref}})}. \end{equation} (18)

For the practical calculation of spectrum functions, $\{ \sigma _{j}^{\text{aux}} \}$ values are selected to be consistent with the real part of $\{ \sigma _{j} \}$ in shifted linear systems, i.e., $\sigma _{j}^{\text{aux}} = \text{Re}\{ \sigma _{j} \}$ (j = 1, …, M). The pseudo-code of the RSCG-SS method is given in Table 1.

Table 1 Pseudo-code of the RSCG-SS method for the evaluation of g(σ_j) = b^†(σ_jI + A)⁻¹b with a Hermitian matrix A.

3. Results and Discussions

TM-L_2,3 XAS spectra of simple TM oxides were computed for testing the performance of the proposed algorithm. MnO and SrTiO₃ were selected as benchmarking systems, in which the TM ions are located at the octahedral (O_h) site. The calculations were made in real space with the cluster models composed of single TM ion and six neighboring oxygen ions, which are embedded into the Madelung potential to mimic the solid-state effects. The atomic positions of these compounds were obtained from the experimental crystalline structures.^34,35) As the first step, fully-relativistic molecular orbital (MO) calculations were performed by solving the four-component Dirac equation within the framework of local density approximation. Relativistic MOs were expressed as linear combinations of numerically generated four-component atomic orbitals (1s-4p for TM and 1s-2p for O).

The ab-initio CI calculations were made using the obtained MOs as basis functions to expand the Hilbert space. In the terminology of quantum chemistry, this method is called as the restricted active space CI (RASCI) calculations. The MOs mainly composed of TM-2p and TM-3d atomic orbitals were selected as active orbitals. The many-electron wavefunctions corresponding to the initial and final state of electric transitions at TM L_2,3-edges were expressed as linear combinations of Slater determinants corresponding to (2p)⁶(3d)ⁿ and (2p)⁵(3d)ⁿ⁺¹ configurations, respectively, where n is the number of TM-3d electrons. The configurations having two or more holes on TM-2p MOs were not considered since the multiplet energies of such configurations are much higher and do not interact with the two configurations described above.

In the case of SrTiO₃, the ligand-to-metal charge transfer (CT) is known to affect the spectral shape of Ti-L_2,3 XAS. Thus, the MOs mainly composed of ligand O-2p atomic orbitals were also treated as active orbitals to consider the CT. Three different calculations were performed in which the maximum number of ligand holes on O-2p orbitals were set to one, two, and three, respectively. Hereafter we refer to these three types of calculations as CT1, CT2, and CT3, respectively. The additional electronic configurations considered in the CT3, for example, are (2p)⁶(O-2p)³⁵(3d)¹, (2p)⁶(O-2p)³⁴(3d)², and (2p)⁶(O-2p)³³(3d)³, for the initial states, and (2p)⁵(O-2p)³⁵(3d)², (2p)⁵(O-2p)³⁴(3d)³, and (2p)⁵(O-2p)³³(3d)⁴ for the final states.

Many-electron Hamiltonian matrices were constructed in terms of Slater determinants corresponding to the electronic configurations described above. Each Hamiltonian matrix was generated in block diagonal form by considering the irreducible representation of C_4h double group, which is the largest abelian subgroup of O_h. Then, a few eigenvalues and eigenvectors in each block matrix were obtained by using the generalized Davidson method. The maximum dimension of block Hamiltonian matrices were 316 for the Mn-L_2,3 XAS of MnO, 1,230 (CT1), 56,310 (CT2), and 270,510 (CT3) for the Ti-L_2,3 XAS of SrTiO₃, respectively. Finally, the PACS were calculated using the RSCG-SS method described in Sec. 2 within the electric dipole approximation.

Figure 1(a) shows the imaginary part of the Green’s function corresponding to the Mn L_2,3-edges of MnO obtained by the RSCG-SS method as a function of incident photon energy ħω and smearing width η. The prescribed accuracy for the RSCG-SS was set to ϵ = 10⁻⁴. The RSCG-SS step converged at the whole energy region displayed in Fig. 1(a) without any numerical instabilities. Figure 1(b) shows the theoretical XAS spectrum calculated by eq. (2) from the cross-section of Im G_A(ω, η) at η = 0.2 eV (dots). The theoretical spectrum obtained by the Fermi’s golden rule (eq. (1)), using the eigenstates obtained by fully-diagonalizing the Hamiltonian matrix, is also shown in Fig. 1(b) (solid line), where the δ-function in eq. (1) is replaced by a Lorentz function, and its full-width at half maximum (FWHM) is 2η = 0.4 eV. Two theoretical spectra obtained by the different algorithms have identical shapes. The results indicate that the RSCG-SS method was implemented successfully.

Fig. 1

Theoretical Mn-L_2,3 XAS of MnO. (a) Imaginary part of Green’s function, Im G_A(ω, η). (b) Cross section of Im G_A(ω, η) at η = 0.2 eV (dots), which is compared with the theoretical spectrum obtained by the full diagonalization method with the same broadening width (solid line).

The seed switching plays an important role in increasing the numerical stability during the RSCG iterations, as mentioned in Sec. 1. To demonstrate this point, we also performed calculations without seed switching. Figure 2 shows the Mn-L_2,3 XAS spectra by the RSCG method with two different reference energies, σ^ref, for the seed system. As can be seen, the spectrum function converged at the whole energy region when taking σ^ref = 647 eV, while it diverged between 640 eV to 655 eV when taking σ^ref = 500 eV after 108 RSCG iterations. We monitored the internal parameters in each iteration and found that the absolute value of the scalar factor, ρ_k(σ_j), which is used for estimating the residual norms in shifted linear systems, became too large to be expressed in the double-precision floating-point number. Such behaviors were observed when reference energy σ^ref was selected far from the absorption edge energy or setting prescribed accuracy ϵ too small. This kind of numerical instabilities was completely suppressed when the seed switching was turned on. In other words, the seed switching acts as the safety mechanism to prevent the numerical overflow of ρ_k(σ_j).

Fig. 2

Theoretical Mn-L_2,3 XAS of MnO calculated by the RSCG method (no seed switching) with two different reference energies, σ^ref = 647 eV, and 500 eV. For the latter case, the Green’s function G_A(ω, η) is diverged after 108 RSCG iterations.

As a test case of applying the proposed algorithm to large Hamiltonian matrices, the Ti-L_2,3 XAS spectra of SrTiO₃ were calculated, including ligand-to-metal CT. Figure 3 shows the theoretical spectra with the different number of ligand holes, namely, CT1, CT2 and CT3, obtained by the RSCG-SS method. The spectra are also compared to the experimental spectra obtained from the literature.³⁶⁾ On these three calculations, the Green’s functions were converged within the desired accuracy. The results clearly demonstrate that the proposed algorithm is stable even in the systems with a sizable Hamiltonian matrix of O(10⁵) dimensions. Comparing three theoretical spectra, one can find that the deviation of spectral shapes between CT1 and CT2: this means that the electronic configuration with two ligand holes contributes to modify the transition energy and intensities of satellite peaks. The results also indicate that CT configurations with more than one ligand holes should be considered for calculations of TM-L_2,3 XAS from high-valence TM compounds. One can also find that both the CT2 and CT3 yield identical spectral shapes, which means that the electronic configuration with three ligand holes has insignificant effects on the multiplet structure and spectral shape.

Fig. 3

Theoretical Ti-L_2,3 XAS of SrTiO₃ obtained by the RSCG-SS method, including the single, double, and triple ligand-to-metal charge transfer, which is referred as CT1, CT2, and CT3, respectively (upper panel). The maximum dimension of block Hamiltonian matrix for these calculations are 1,230 for CT1, 56,310 for CT2, and 270,510 for CT3, respectively. These are compared with the experimental spectrum from the literature³⁶⁾ (lower panel).

Finally, we would mention about the computational efficiency of the RSCG-SS method. In each iteration, a single matrix-vector multiplication operation is required to compute $A\boldsymbol{{r}}_{k}^{\text{ref}}$, which appears in eqs. (5) and (7). Other operations for updating the residual vector $\boldsymbol{{r}}_{k}^{\text{ref}}$ are done in O(N), and for evaluation of Green’s function g(σ₁), …, g(σ_M) requires O(M) in each iteration. The total computational time of the RSCG-SS method scales as n_iter (K + O(N) + O(M)), where K is the computational complexity of the matrix-vector multiplication, and n_iter is the number of iterations until convergence. Both K and n_iter depends on the character of matrix A. If n_iterK is much smaller than O(N³), the computational time of the RSCG-SS method would be faster than that of full-diagonalization method. In the case that a matrix A is sparse (this means, most of the matrix elements are zero), or it has a special structure such that the matrix-vector product can be performed efficiently, the RSCG-SS algorithm can be efficient.

In the RSCG-SS method, matrix A is required only for computing the matrix-vector multiplication. At this step, only the non-zero matrix elements are required. The memory space required to store matrix elements can be significantly reduced if matrix A is sparse. Alternatively, we can compute matrix elements of A on-demand in each iteration to further reduce the memory consumption with the expense of computational time. This enables us to apply the RSCG-SS method for a huge matrix whose dimension is of O(10⁸) or more, which never be handled by the full-diagonalization method.

In the CI calculations used in this work, we can utilize the sparseness of a Hamiltonian matrix, which makes the RSCG-SS method efficient. Let us show some benchmarking results of CT1, CT2, and CT3 calculations for SrTiO₃. All the calculations were performed on a laptop computer with Intel^® Core ™ i7 processor @ 3.3 GHz. All the floating-point arithmetic was done in double-precision. As the CI Hamiltonian matrix is sparse, matrix elements were stored in the compressed sparse row (CSR) format. The computational time was about 3 minutes for CT1, 20 minutes for CT2, and two hours for CT3, respectively. The memory space used in these calculations was about 1 MB for CT1, 160 MB for CT2, and 1 GB for CT3, respectively. Much larger memory space is required in the case full-diagonalization approach because whole matrix elements have to be stored on memory for dense matrix decomposition. The estimated memory required is about 20 MB for CT1, 50 GB for CT2, and 1 TB for CT3, respectively. The full-diagonalization method is unfeasible for the CT3, even using the modern workstations. These results clearly show that the RSCG-SS method is quite efficient for calculating Green’s functions corresponding to a large Hamiltonian matrix. We expect the RSCG-SS method can be applied to a larger matrix of O(10⁶) dimensions.

4. Conclusion

In this study, a new iterative algorithm, called the RSCG-SS method, for evaluating spectrum function in eq. (3) was developed, in which the automatic seed switching mechanism was introduced into the RSCG method. The spectrum functions with different energies, ħω, and smearing factors, η, can be obtained simultaneously by using the proposed method. The RSCG-SS method was applied to calculate TM-L_2,3 XAS spectra of simple TM oxides as benchmarks. The theoretical spectra obtained by the RSCG-SS method were identical to that obtained by the full-diagonalization method in accordance with the Fermi’s golden rule. The RSCG-SS method was shown to be numerically stable even in systems of a sizable Hamiltonian matrix of O(10⁵) dimension. The RSCG-SS is suitable for simulating the spectra from larger systems, such as metal clusters and TM complexes whose ligands are composed of multiple of atoms. Further research in this direction is underway.

Acknowledgments

This work was supported by PRESTO, Grant no. JPMJPR16N1 16815006 from Japan Science and Technology Agency (JST), and the Grants-in-Aid for Scientific Research, Grant no. 19H00818 from the Japan Society for the Promotion of Science (JSPS), the Ministry of Education, Culture, Sports, Science and Technology of Japan.

REFERENCES

• 1)  G. Bunker: Introduction to XAFS: A Practical Guide to X-Ray Absorption Fine Structure Spectroscopy, 1st ed., (Cambridge University Press, Cambridge, 2010).
• 2)  S. Calvin: XAFS for Everyone, 1st ed., (Taylor & Francis, Boca Raton, 2013).
• 3)   W.-Y.  Ching and  P.  Rulis: J. Phys. Condens. Matter 21 (2009) 104202.
• 4)   S.-D.  Mo and  W.Y.  Ching: Phys. Rev. B 62 (2000) 7901–7907.
• 5)   C.  Hébert: Micron 38 (2007) 12–28.
• 6)   T.  Mizoguchi,  I.  Tanaka,  S.  Yoshioka,  M.  Kunisu,  T.  Yamamoto and  W.Y.  Ching: Phys. Rev. B 70 (2004) 045103.
• 7)   T.  Mizoguchi,  I.  Tanaka,  S.-P.  Gao and  C.J.  Pickard: J. Phys. Condens. Matter 21 (2009) 104204.
• 8)   S.-P.  Gao,  C.J.  Pickard,  A.  Perlov and  V.  Milman: J. Phys. Condens. Matter 21 (2009) 104203.
• 9)   J.J.  Rehr: Radiat. Phys. Chem. 75 (2006) 1547–1558.
• 10)   J.J.  Rehr and  R.C.  Albers: Rev. Mod. Phys. 72 (2000) 621–654.
• 11)   A.  Filipponi and  A.  Di Cicco: Phys. Rev. B 52 (1995) 15135–15149.
• 12)   Y.  Joly: Phys. Rev. B 63 (2001) 125120.
• 13)   E.L.  Shirley: J. Electron Spectrosc. Relat. Phenom. 144–147 (2005) 1187–1190.
• 14)   R.  Laskowski and  P.  Blaha: Phys. Rev. B 82 (2010) 205104.
• 15)   W.  Olovsson,  I.  Tanaka,  T.  Mizoguchi,  P.  Puschnig and  C.  Ambrosch-Draxl: Phys. Rev. B 79 (2009) 041102.
• 16)   K.  Tomita,  T.  Miyata,  W.  Olovsson and  T.  Mizoguchi: J. Phys. Chem. C 120 (2016) 9036–9042.
• 17)  F. de Groot and A. Kotani: Core Level Spectroscopy of Solids, (CRC Press, Boca Raton, 2008).
• 18)   F.M.F.  de Groot: J. Electron Spectrosc. Relat. Phenom. 67 (1994) 529–622.
• 19)   H.  Ikeno,  F.M.F.  de Groot,  E.  Stavitski and  I.  Tanaka: J. Phys. Condens. Matter 21 (2009) 104208.
• 20)   H.  Ikeno,  T.  Mizoguchi and  I.  Tanaka: Phys. Rev. B 83 (2011) 155107.
• 21)   H.  Ikeno: Mater. Trans. 56 (2015) 1448–1451.
• 22)   R.V.  Pinjari,  M.G.  Delcey,  M.  Guo,  M.  Odelius and  M.  Lundberg: J. Chem. Phys. 141 (2014) 124116.
• 23)   M.W.  Haverkort,  M.  Zwierzycki and  O.K.  Andersen: Phys. Rev. B 85 (2012) 165113.
• 24)   M.W.  Haverkort: J. Phys. Conf. Ser. 712 (2016) 012001.
• 25)   H.  Ikeno and  T.  Mizoguchi: Microscopy 66(5) (2017) 305–327.
• 26)   T.  Mizoguchi,  W.  Olovsson,  H.  Ikeno and  I.  Tanaka: Micron 41 (2010) 695–709.
• 27)   I.  Tanaka,  T.  Mizoguchi and  T.  Yamamoto: J. Am. Ceram. Soc. 88 (2005) 2013–2029.
• 28)   M.  Urasaki and  H.  Ikeno: Mater. Trans. 59 (2018) 1057–1061.
• 29)   J.J.  Rehr,  J.J.  Kas,  M.P.  Prange,  A.P.  Sorini,  Y.  Takimoto and  F.  Vila: C. R. Phys. 10 (2009) 548–559.
• 30)   P.  Giannozzi et al. J. Phys. Condens. Matter 21 (2009) 395502.
• 31)   Y.  Nagai,  Y.  Shinohara,  Y.  Futamura and  T.  Sakurai: J. Phys. Soc. Jpn. 86 (2017) 014708.
• 32)   S.  Yamamoto,  T.  Sogabe,  T.  Hoshino,  S.-L.  Zhang and  T.  Fujiwara: J. Phys. Soc. Jpn. 77 (2008) 114713.
• 33)  J. Liesen and Z. Strakos: Krylov Subspace Methods: Principles and Analysis, (Oxford University Press, Oxford, 2015).
• 34)   A.N.  Christensen: Acta Chem. Scand. 44 (1990) 851–852.
• 35)   M.  Kuriyama and  S.  Hosoya: J. Phys. Soc. Jpn. 17 (1962) 1022–1029.
• 36)   Y.  Uehara,  D.W.  Lindle,  T.A.  Callcott,  L.T.  Terminello,  F.J.  Himpsel,  D.L.  Ederer,  J.H.  Underwood,  E.M.  Gullikson and  R.C.C.  Perera: Appl. Phys., A Mater. Sci. Process. 65 (1997) 179–182.