EDRIXS: An open source toolkit for simulating spectra of resonant inelastic x-ray scattering

TL;DR

EDRIXS is an open source software toolkit that simplifies the simulation of resonant inelastic x-ray scattering spectra in strongly correlated materials, aiding researchers in understanding complex experimental data.

Contribution

The paper introduces EDRIXS, a comprehensive open source toolkit for simulating RIXS spectra using exact diagonalization and modern computational techniques.

Findings

Efficient simulation of RIXS spectra for various models

Integration with DFT, Wannier90, and DMFT data

Demonstrated usage through three example cases

Abstract

Resonant inelastic x-ray scattering (RIXS) has become a very powerful experimental technique to probe a broad range of intrinsic elementary excitations, for example, from low energy phonons and (bi-)magnons to high energy $d$-$d$, charge-transfer and plasmons excitations in strongly correlated electronic systems. Due to the complexity of the RIXS cross-section and strong core-hole effects, theoretical simulation of the experimental RIXS spectra is still a difficult task which hampers the understanding of RIXS spectra and the development of the RIXS technique. In this paper, we present an open source toolkit (dubbed EDRIXS) to facilitate the simulations of RIXS spectra of strongly correlated materials based on exact diagonalization (ED) of certain model Hamiltonians. The model Hamiltonian can be from a single atom, small cluster or Anderson impurity model, with model parameters from density functional theory plus Wannier90 or dynamical mean-field theory calculations. The spectra of x-ray absorption spectroscopy (XAS) and RIXS are then calculated using Krylov subspace techniques. This toolkit contains highly efficient ED, XAS and RIXS solvers written in modern Fortran 90 language and a convenient Python library used to prepare inputs and set up calculations. We first give a short introduction to RIXS spectroscopy, and then we discuss the implementation details of this toolkit. Finally, we show three examples to demonstrate its usage.

Figures

Tables

Table 1. Table 1: Requirement of tools and external libraries for building EDRIXS code.

Tools and Libraries	Recommendation
Fortran compiler	Intel’s ifort or GNU’s gfortran
MPI	OpenMPI or MPICH
BLAS/LAPACK	MKL or OpenBLAS [49]
Parallel ARPACK	arpack-ng [50], v3.6.2 or higher
Python3	v3.6 or higher
Numpy	v1.15 or higher
Scipy	v1.1 or higher
Sympy	v1.3 or higher
Matplotlib	v3.0 or higher
Sphinx	v1.7 or higher
Numpydoc	v0.8 or higher

Table 2. Table 2: Input and output files for ed.x .

	File Name	Description
Inputs	config.in	set up control parameters
	fock_i.in	Fock basis $|I⟩$ for ${\hat{H}}_i$
	hopping_i.in	$t_{\alpha ,\beta }$ terms in ${\hat{H}}_i$
	coulomb_i.in	$U_{\alpha ,\beta ,\gamma ,\delta }$ terms in ${\hat{H}}_i$
Outputs	standard outputs	log file
	eigvals.dat	eigenvalues $E_i$
	eigvec.$i$	the $i$-th ground state $|{\Gamma }_i⟩$
	denmat.dat	density matrix: $⟨{\Gamma }_i\left|{\hat{f}}_{\alpha }^{†}{\hat{f}}_{\beta }\right|{\Gamma }_i⟩$

Table 3. Table 3: Input and output files for xas.x .

	File Name	Description
Inputs	config.in	set up control parameters
	fock_i.in	Fock basis $|I⟩$ for ${\hat{H}}_i$
	fock_n.in	Fock basis $|I⟩$ for ${\hat{H}}_n$
	hopping_n.in	$t_{\alpha ,\beta }$ terms in ${\hat{H}}_n$
	coulomb_n.in	$U_{\alpha ,\beta ,\gamma ,\delta }$ terms in ${\hat{H}}_n$
	transop_xas.in	transition operator ${\widehat{𝒟}}_i$
	eigvec.$i$	the $i$-th ground state $|{\Gamma }_i⟩$ of ${\hat{H}}_i$
Outputs	standard outputs	log file
	xas_poles.$i$	XAS data for the $i$-th ground state

Table 4. Table 4: Input and output files for rixs.x .

	File Name	Description
Inputs	config.in	set up control parameters
	fock_i.in	Fock basis $|I⟩$ for ${\hat{H}}_i$
	fock_n.in	Fock basis $|I⟩$ for ${\hat{H}}_n$
	fock_f.in	Fock basis $|I⟩$ for ${\hat{H}}_i$
	hopping_i.in	$t_{\alpha ,\beta }$ terms in ${\hat{H}}_i$
	hopping_n.in	$t_{\alpha ,\beta }$ terms in ${\hat{H}}_n$
	coulomb_i.in	$U_{\alpha ,\beta ,\gamma ,\delta }$ terms in ${\hat{H}}_i$
	coulomb_n.in	$U_{\alpha ,\beta ,\gamma ,\delta }$ terms in ${\hat{H}}_n$
	transop_rixs_i.in	transition operator ${\widehat{𝒟}}_i$
	transop_rixs_f.in	transition operator ${\widehat{𝒟}}_f^{†}$
	eigvec.$i$	the $i$-th ground state $|{\Gamma }_i⟩$ of ${\hat{H}}_i$
Outputs	standard outputs	log file
	rixs_poles.$i$	RIXS data for the $i$-th ground state

Equations

$$I_{\text{XAS}}({\omega }_{\text{in}},{\overrightarrow{ϵ}}_i)=\sum _i\frac{1}{Z}{e}^{-\frac{E_i}{K_{B}T}}\sum _n{\left|⟨n⟩{\hat{\mathcal{D}}}_{i}i\right|}^2\frac{{\Gamma }_c/\pi }{({\omega }_{\text{in}}-E_n+E_i{)}^2+{\Gamma }_c^2},$$

$$I_{\text{RIXS}}({\omega }_{\text{in}},\omega ,{\vec{k}}_i,{\vec{k}}_f,{\overrightarrow{ϵ}}_i,{\overrightarrow{ϵ}}_f)$$

$${\hat{\mathcal{D}}}_i$$

$${\hat{\mathcal{D}}}_f^{†}$$

$${\hat{T}}_x^i$$

$${\hat{T}}_x^f$$

$${\hat{H}}_i=\sum _{\alpha ,\beta }{t}_{\alpha ,\beta }{\hat{f}}_{\alpha }^{†}{\hat{f}}_{\beta }+\sum _{\alpha ,\beta ,\gamma ,\delta }{U}_{\alpha ,\beta ,\gamma ,\delta }{\hat{f}}_{\alpha }^{†}{\hat{f}}_{\beta }^{†}{\hat{f}}_{\gamma }{\hat{f}}_{\delta },$$

$${\hat{H}}_n={\hat{H}}_i+{\hat{V}}_{\text{core-hole}}+{\hat{H}}_{\text{core}},$$

$${\hat{H}}_i={\hat{H}}_{\text{imp}}+\sum _{\alpha ,l}{E}_{\alpha ,l}{\hat{c}}_{\alpha ,l}^{†}{\hat{c}}_{\alpha ,l}+\sum _{\alpha ,l}{V}_{\alpha }^l{\hat{f}}_{\alpha }^{†}{\hat{c}}_{\alpha ,l}+h.c.,$$

$${\Delta }_{\alpha }(i{\omega }_n)=\sum _{l=1}^N\frac{{\left|V_{\alpha }^l\right|}^2}{i{\omega }_n-E_{\alpha ,l}}.$$

$$\left({\omega }_{\text{in}}-{\hat{H}}_n+E_i+i{\Gamma }_c\right)|x⟩=|b⟩,$$

$$I_{\text{RIXS}}({\omega }_{\text{in}},\omega )=⟨F⟩\frac{1}{\omega -{\hat{H}}_i+E_i+i\Gamma }F.$$

$$\overset{1-6}{\overbrace{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}110100}}}\;\overset{7-12}{\overbrace{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}111011}}}$$

$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\times 2^{1-1}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\times 2^{2-1}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\times 2^{3-1}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}\times 2^{4-1}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\times 2^{5-1}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}\times 2^{6-1}=11.$$

$$\hat{O}=\sum _{\alpha ,\beta }{t}_{\alpha ,\beta }{\hat{f}}_{\alpha }^{†}{\hat{f}}_{\beta }+\sum _{\alpha ,\beta ,\gamma ,\delta }{U}_{\alpha ,\beta ,\gamma ,\delta }{\hat{f}}_{\alpha }^{†}{\hat{f}}_{\beta }^{†}{\hat{f}}_{\gamma }{\hat{f}}_{\delta },$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program title: EDRIXS

Catalogue identifier: TO BE DONE

Program summary URL: TO BE DONE

Program obtainable from: TO BE DONE

Licensing provisions: GNU General Public Licence 3.0

No. of lines in distributed program, including test data, etc.: 28,365 lines

No. of bytes in distributed program, including test data, etc.: 412,888 bytes

Distribution format: tar.gz

Programming language: Fortran 90 and Python3

Computer: Desktop PC, laptop, high performance computing cluster

Operating system: Unix, Linux, Mac OS X

Has the code been vectorised or parallelized?: Yes, it is parallelized by MPI

RAM: Heavily depends on the complexity of the problem

Classification: 7.3

External routines/libraries used: BLAS, LAPACK, Parallel ARPACK, numpy, scipy, sympy, matplotlib, sphinx, numpydoc

Nature of problem: Simulating the spectra of resonant inelastic x-ray scattering for strongly correlated electronic systems.

Solution method: Exact diagonalization of model Hamiltonians and Krylov subspace methods.

Running time: Heavily depends on the complexity of the problem.

1 Introduction

Directly measuring the elementary excitations in materials through experimental spectroscopy is a fundamental task in condensed matter physics. For example, angle-resolved photoemission spectroscopy (ARPES) [1] is a very powerful technique to measure the single particle spectrum function $A(\vec{k},\omega)$ (bands) in both weakly and strongly correlated materials. In recent years, resonant inelastic x-ray scattering (RIXS) has also emerged as a key probe of the elementary excitations in materials, especially for strongly correlated systems [2]. RIXS is a photon-in photon-out spectroscopy. The energy of incident photon is tuned at a specific resonant edge to excite a deep core electron to the valence shell and create a core-hole. After a very short time, the core-hole is refilled by a valence electron and photon is emitted. By measuring the change in energy, momentum and polarization of the scattered photon, one can obtain detailed information on the nature of the underlying excitations. Compared to other spectroscopies, RIXS has various unique features and advantages [2]. First, it can measure a broad range of elementary excitations, for example, from low energy phonons and (bi-)magnons to high energy local $d$-$d$, charge transfer (CT) and plasmon excitations. Second, since x-rays carry appreciable momentum, the dispersion of low-energy excitations in solids can be probed. Third, there are fewer restrictions to sample surface quality or volume compared to probes such as ARPES and neutron scattering, which enables, for instance, studies of thin films and very small single crystals. Finally, it is element and orbital specific, bulk sensitive and polarization dependent. In recent years, great progress in the technical development of RIXS has occurred due to high-brilliance synchrotron light sources and advanced photon detection instrumentation, leading to its wide use to study elementary excitations in strongly correlated materials [2, 3, 4, 5, 6].

Given the RIXS sensitivity to a wide range of elementary excitations, theoretical simulation of the RIXS cross-section is often very relevant in the interpretation of the experimental RIXS spectra and its underlying physics. The RIXS cross-section is described by the Kramers-Heisenberg formula, which is a result of second-order perturbation theory. The RIXS cross-section can be interpreted as a four-particle correlation function involving both valence and core electrons [7, 8] and can be in principle determined by some first-principles numerical methods. Nowadays, density functional theory plus dynamical mean-field theory (DFT+DMFT) [9, 10] with continuous-time quantum Monte-Carlo (CTQMC) [11] as the impurity solver is the most powerful and accurate method for first-principle calculations of physical properites, such as the single particle spectrum function, of strongly correlated materials, so it is a desirable method to compute RIXS spectra. Unfortunately, it is still a very challenging mission to determine the RIXS cross-section by this method with the current available numerical techniques due to the complexity of this correlation function and the strong unknown core-hole potential. The difficulty in performing theoretical RIXS simulations hampers the progress of the understanding of RIXS spectra and the development of the RIXS technique. Nowadays, most of RIXS simulations are based on the exact diagonalization (ED) of Hamiltonians of small size models, that capture sort of the properites of the real materials, such as single atom, small cluster and Anderson impurity model (AIM) [12]. Currently, the most widely used RIXS simulation codes based on ED algorithms include CTM4RIXS [13] and Quanty [14]. The former is based on an old atomic code called Cowan’s code [15] with very limited functionalities to perform XAS and RIXS calculations. The latter is based on a lightweight programming language called Lua [16] and has much richer functionalities. However, the core of the Quanty code is not open source, so it cannot be freely modified by users to implement functionalities that are not available. Many RIXS simulations based on ED [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38] or DMRG [37, 39, 40] have been also reported. In Refs. [25, 28, 30], the ED code is based on the parallel ARPACK library. In Refs. [36], a powerful ED solver based on quantum chemistry algorithm is used to diagonalize the AIM from a converged DFT+DMFT calculations with many bath sites and continuous charge excitations have been found in high-valence transition-metal oxides in their simulated RIXS spectra. However, currently these codes are not open source.

In this work, we present the EDRIXS open source toolkit that is aimed at facilitating RIXS simulations based on ED algorithms. EDRIXS is designed as a post-processing tool for the open source COMSCOPE project [41, 42]. It implements parallel ED solvers with very high efficiency based on the standard Lanczos algorithms and parallel version of ARPACK library. The XAS and RIXS spectra are then calculated parallelly using Krylov subspace algorithms. The core components of this toolkit are implemented using the modern Fortran 90 language and the parallelism is achieved by MPI. The code has a layered structure, thus new functionalities can be easily implemented by users who want to perform other types of x-ray scattering simulations. An application programming interface (API) is implemented by the popular Python3 language to prepare inputs and set up calculations. A pure Python ED solver is also implemented in the API for small size problem. The rest of this paper is organized as follows: In section 2, we introduce the basic theory and methods of RIXS simulations. In section 3, we explain the implementation of the code in detail. In section 4, we show how to install and use this code. In section 5, we show three examples to demonstrate the usage of this code. In section 6, we discuss the plans of future development of EDRIXS.

2 Basic theory and methods of RIXS simulation

2.1 Geometry of RIXS experiment

A typical geometry of RIXS experiment is illustrated in Fig. 1. $\vec{k}_{i}$ ($\vec{\epsilon}_{i}$) and $\vec{k}_{f}$ ($\vec{\epsilon}_{f}$) are the wave (polarization) vectors of the incident and scattered x-ray, respectively. $\omega_{\text{in}}$ and $\omega_{\text{out}}$ are the energy of incident and scattered x-ray, respectively. $\alpha$ and $\beta$ are the angles between the polarization vectors and the scattering plane. We call it $\pi$-polarization when the polarization vector is parallel to the scattering plane ($\alpha,\beta=0$), and $\sigma$-polarization when the polarization vector is perpendicular to the scattering plane ($\alpha,\beta=\pi/2$). $\theta_{\text{in}}$ and $\theta_{\text{out}}$ are the incident and scattered angles, respectively. $\phi$ is the azimuthal angle defined with respect to $x$-axis.

2.2 XAS and RIXS cross-section

The XAS cross-section is described by,

[TABLE]

and the RIXS cross-section is described by the Kramers-Heisenberg formula [2],

[TABLE]

where, $\Ket{i}$ are the ground states of the Hamiltonian $\hat{H}_{i}$ describing the valence electrons. $T$ is temperature and $K_{B}$ is the Boltzmann factor. $Z=\sum_{i}e^{-\frac{E_{i}}{K_{B}T}}$ is the partition function. $\Ket{f}$ are the excited eigenstates of $\hat{H}_{i}$. $\Ket{n}$ are the eigenstates of intermediate Hamiltonian $\hat{H}_{n}$ which includes a deep core-hole and one more electron in the valence shell after the x-ray absorption process. $\omega_{\text{in}}$ is the energy of the incident x-ray which is tuned to be at a specific resonant edge. $\omega$ is the x-ray energy loss, i.e. the energy difference between the incident and scattered x-ray. $\Gamma_{c}$ and $\Gamma$ are the lifetime broadening of the intermediate (with a core hole) and final states, respectively.

$\hat{\mathcal{D}}_{i}$ and $\hat{\mathcal{D}}^{\dagger}_{f}$ are the transition operators for x-ray absorption and emission process, respectively. They can be generally written as,

[TABLE]

where, $P_{a}^{i}$ and $P_{a}^{f}$ involve the polarization and wave vector of x-ray. For dipolar transition, they are: $\epsilon_{x}^{i},\epsilon_{y}^{i},\epsilon_{z}^{i}$ and $\epsilon_{x}^{f},\epsilon_{y}^{f},\epsilon_{z}^{f}$. $\hat{T}_{a}^{i}$ and $\hat{T}_{a}^{f}$ are components of transition operators. For dipolar transition, they are:

[TABLE]

where, $\vec{R}$ is site-index, and $\hat{x}_{R}$, $\hat{y}_{R}$ and $\hat{z}_{R}$ are position operators of electrons bound to site-$R$. These are measured with respect to the center of site-$R$.

2.3 Model Hamiltonian

$\hat{H}_{i}$ is a model Hamiltonian for valence electrons from single atom, cluster or Anderson impurity models. It can be generally witten as,

[TABLE]

where, $t_{\alpha,\beta}$ usually contain crystal field (CF) splitting terms, spin-orbit coupling (SOC) terms, hopping terms between different sites. The CF and hopping terms can be obtained from a first-principle DFT+Wannier90 calculations. SOC can be obtained from DFT calculations or single atomic calculation, for example, from Cowan’s atomic code [15]. $U_{\alpha,\beta,\gamma,\delta}$ are the rank-4 tensors of the on-site Coulomb interaction and are parameterized by Slater integrals, which can be obtained from an atomic calculations.

$\hat{H}_{n}$ is the intermediate Hamiltonian with a core-hole, which can be written as,

[TABLE]

where, the core-hole potential $\hat{V}_{\text{core-hole}}$ is simulated at atomic level by adding Coulomb interaction between valence electrons and core-hole. $\hat{H}_{\text{core}}$ is the Hamiltonian of the core electrons.

To simulate RIXS spectra based on DFT+DMFT calculations, we first perform a DFT+DMFT calculation to get converged hybridization function $\Delta_{\alpha}(i\omega_{n})$ and then construct an Anderson impurity model $\hat{H}_{i}$ from it, which can be written as,

[TABLE]

where, $\hat{H}_{\text{imp}}$ is the Hamiltonian for the impurity site, which is taken from the atomic problem in DFT+DMFT calculation. $E_{\alpha,l}$ is the energy level of the $l$-th bath site with orbital $\alpha$, $V^{l}_{\alpha}$ are the hybridization strength between localized impurity electrons ($\hat{f}^{\alpha}$) and conducted bath electrons ($\hat{c}_{\alpha,l}$). $E_{\alpha,l}$ and $V_{\alpha}^{l}$ are obtained by fitting the hybridization function,

[TABLE]

The core-hole potential is only added to the impurity site in $\hat{H}_{n}$.

2.4 RIXS simulation algorithms

We use similar RIXS simulation algorithms as that used in Ref. [25, 36]. The flow diagram of the algorithm to calculate the RIXS cross-section is illustrated in Fig. 2 and briefly described as follows:

Use an ED solver to diagonalize $\hat{H}_{i}$ to get the ground states or a few low-energy excited states $\Ket{i}(E_{i})$. Both standard Lanczos algorithm [43] without re-orthogonalization and implicitly restarted Arnoldi algorithm [44, 45] are used in ED solvers.

• 2.

Apply the transition operator for the x-ray absorption process on $\Ket{i}$ to get $\Ket{b}=\hat{\mathcal{D}}_{i}\Ket{i}$.

• 3.

Solve the following large-scale sparse symmetric linear equations,

[TABLE]

by minimum residual (MINRES) method [46].

• 4.

Apply the transition operator for the x-ray emission process on $\Ket{x}$ to get $\Ket{F}=\hat{\mathcal{D}}^{\dagger}_{f}\Ket{x}$.

• 5.

Build Krylov space of $\hat{H}_{i}$ and use continued fraction technique to get the RIXS spectra,

[TABLE]

3 Implementations and optimizations

In this section, we describe the details of implementation and optimization of EDRIXS.

3.1 Development platform

The key components of the EDRIXS code are developed with the modern Fortran 90 language. Intel’s ifort and GNU’s gfortran compilers have been tested. The code is parallelized by MPI. An interface that is used for pre- and post-processing is written by Python3 language. Git is used as the version control system.

3.2 Fock basis

We implement EDRIXS code based on the second-quantization language. We first define a single particle basis $\alpha$, then define Fock basis $\Ket{I}$ with respect to it. In the code, all the orbitals of valence electrons are gathered together and put in front of all the orbitals of core electrons (i.e. indices of valence orbitals are always smaller than those of core orbitals). A Fock state can be represented by a binary number with digital 1 or 0, where 1 means that the orbital is occupied and 0 means that the orbital is empty. For example,

[TABLE]

where, the red parts with orbital indices from 1 to 6 are for valence electrons and the blue parts with orbital indices from 7 to 12 are for core electrons. In this example, valence orbitals 1, 2 and 4 are occupied, and one core orbital 10 is empty. The above binary format is used for the ED solver implemented in Python API. For ed.x, xas.x and rixs.x, we need to input the Fock basis with their corresponding decimal formats, but only for the parts of valence orbitals. The core orbitals will be explicitly considered within these three programs. For the above example (Eqn. 13), the corresponding decimal number is computed as following,

[TABLE]

3.3 Format of sparse matrices

In second-quantization form, a general operator can be written as,

[TABLE]

where, we call $t_{\alpha,\beta}$ hopping terms and $U_{\alpha,\beta,\gamma,\delta}$ rank-4 interaction terms. The matrix elements of $\hat{O}$ in the Fock basis $\Ket{I}$ are $\Braket{I^{\prime}}{\hat{O}}{I}$. These matrices are large-scale sparse matrices and the number of nonzero elements is proportional to its dimension $n$, so we store them in compressed sparse row (CSR) format to save memory. To speed up the building of operators, only nonzeros of $t_{\alpha,\beta}$ and $U_{\alpha,\beta,\gamma,\delta}$ are input.

3.4 Structure of the code

The structure of the EDRIXS code can be divided into four layers, which is illustrated in Fig. 3. The key components of EDRIXS code involve two basic operations: building matrices form of operators in Fock basis and performing sparse matrix-vector multiplication (MVM), which are the most time-consuming parts, so we implement subroutines to parallelly do these tasks in the first layer. The rows of a sparse matrix $M$ and elements of a vector $V$ are equally distributed among all the processors. Each processor only needs to build parts of the rows of the operator and do parts of the MVM, and all the processors can work parallelly. The only MPI communication needed for MVM is to send and receive the vector $V$, which is a bottleneck of the code. The speed of this MPI communication heavily depends on the networking of the cluster. Based on the first layer, we implement subroutines to parallelly build tridiagonal Krylov space and solve symmetric linear equations using MINRES methods in the second layer. The key operations of these two subroutines are MVM.

Above the first and second layers, we implement the main functionalities: ed.x, xas.x and rixs.x to do RIXS simulations in the third layer. Three types of ED solvers are implemented. The first one is to use the LAPACK’s subroutine ZHEEV to diagonalize a Hamiltonian to get all the eigenstates, which is useful for small size ($n<10,000$) problem. The second one is to use standard Lanczos algorithm without re-orthogonalization to get only one ground state, which is not very accurate due to the orthogonality problem, especially for degenerate ground states. However, it runs very fast and is suitable for roughly searching the ground state of a Hamiltonian in many different subspaces. The third one is to use the parallel version of the ARPACK library to get few lowest eigenstates, which is very accurate because it uses implicitly restarted Arnoldi algorithm to deal with the orthogonalization issue, but its efficiency is lower than the second one. We will use this solver in real RIXS simulations. Based on the first and second layers, we can easily extend new functionalities to do more calculations, for example, calculating the expected value of an operator in the ground state and other correlation functions.

In the fourth layer, a Python interface is implemented to provide some application programming interfaces (APIs) for preparing inputs for ed.x, xas.x and rixs.x and setting up calculations. Besides, a pure python ED solver based on second-quantization is also implemented for small size problem, so the RIXS simulations can be done easily with pure Python code. Currently, we have not implemented the interface with the Cowan’s code for calculating Slater integrals. One should manually run Cowan’s code and read Slater integrals. We refer users to its official webpage [47] for more information.

4 Installation and usage

In this section, we will show how to install and use EDRIXS.

4.1 Get EDRIXS

EDRIXS is an open source code and is released under the GNU General Public Licence 3.0 (GPL). The readers who are interested in it can write a letter to the authors to request a copy of the newest version of EDRIXS source code, or they can download it directly from the public code repository:

$ git clone https://github.com/shenmidelin/edrixs.git.

where $ is the command line prompt. One can use the command

$ git pull

to get the latest updates of the code.

4.2 Build EDRIXS

To build the source code and documentation of EDRIXS, some tools and external libraries are required, which are listed in Table 1. Most of the libraries can be easily installed using Anaconda [48]. There may be problems when using gfortran with MKL, so we strongly recommend gfortran+OpenBLAS. Note that the arpack-ng library should be compiled with the same Fortran compiler and BLAS/LAPACK libraries.

We assume that one downloads the source code into the directory EDRIXS_DIR. One then go to the edrixs/src/fortran directory and edit the make.sys file to configure the compiling settings.

$ cd edrixs/src/fortran
$ cp make.sys.ifort make.sys  (or cp make.sys.gfortran make.sys)
$ vim make.sys

where, in the make.sys file, one needs to set F90 as

F90 = mpif90

and set LIBS, for example,

LIBS = -L${MKLROOT}/lib/intel64 -lmkl_core -lmkl_sequential -lmkl_rt \\
       -L${ARPACK_DIR}/lib/ -lparpack -larpack

where, we assume the ARPACK library is installed in ARPACK_DIR. Then, we type

$ make
$ make install

to compile and install the executable files ed.x, xas.x and rixs.x into bin directory. Then, one should add the following two lines in .bashrc or .bash_profile file,

export PATH=${EDRIXS_DIR}/edrixs/bin:$PATH
export PYTHONPATH=${EDRIXS_DIR}/edrixs/src/python:$PYTHONPATH

To compile the documentation of the Python APIs, one should

$ cd edrixs/docs
$ mkdir build
$ sphinx-build -b html source build
$ make html

and open the file

${EDRIXS_DIR}/edrixs/docs/build/index.html

in a browser to read the documentation of the Python APIs.

4.3 Use EDRIXS

The typical procedures of the usage of EDRIXS are as follows,

Use Python APIs to set up the parameters of Hamiltonian, transition operators, Fock basis.

• 2.

Run ed.x, xas.x and rixs.x in order. These three Fortran executable files can be launched manually, for example,

mpirun -np 16 ed.x > log.dat

or through the Python subprocess call, for example,

import subprocess
subprocess.check_call([‘‘mpirun", ‘‘-np", ‘‘16", ‘‘ed.x"])

• 3.

Use Python APIs to post-process the results and plot XAS or RIXS spectra.

The required input and output files for ed.x, xas.x and rixs.x are listed in Table 2, Table 3 and Table 4, respectively. More details about the input and output files are described in file user_guide.pdf in docs directory.

5 Examples

In this section, we show three examples to demonstrate the usage of EDRIXS.

5.1 Hello RIXS: multiplets $d$-$d$ excitations in a single atom case Ni ($d^{8}$)

In this example, we use pure Python code to do ED calculation and calculate XAS and RIXS spectra step by step. The main purpose is to show how the Python APIs are used. This example is a resonant x-ray scattering at Ni-$L_{2/3}$ edge ($2p_{1/2,3/2}\rightarrow 3d$ transition). Ni has a $d^{8}$ configuration with a tetragonal CF environment [3]. The Python codes are shown in the Listing 1.

The calculated XAS spectrum and a RIXS map at $L_{3}$ edge are plotted in Fig. 4. More results can be found in Ref. [3].

5.2 Two-site Ir-Ir cluster: dimer excitations

In this section, we show an example of two-sites Ir-Ir cluster model, where two IrO6 octahedras share their face, see Fig. 5. This type of structure can be found in materials, such as Ba5AlIr2O11 [51, 52] and the 6H-hexagonal oxides Ba${}_{3}AB_{2}$O9 ($A$=In, Y, Lu, Na and $B$=Ru, Ir) [53, 54, 55]. Such face sharing structure leads to strong hopping between these two Ir sites, and dimer excitations have been observed in experimental RIXS spectra and confirmed by a two-site model simulation [56]. Here, we show how to do the calculations step by step. We will simulate RIXS spectra at Ir-$L_{3}$ edge for the cases with and without hopping between these two Ir sites based on a $t_{2g}$ two-sites model. More details of the methods and the Hamiltonian of the on-site Coulomb interaction, spin-orbit coupling (SOC), trigonal crystal field and the hoppings between the two Ir sites can be found in Ref. [56].

We first go to the examples directory and run the Python script,

$ cd ${EDRIXS_DIR}/examples/cpc/two_site_cluster
$ ./get_inputs.py

to do the following things,

It first makes four directories: ed for performing ED calculation, xas for performing XAS calculation, rixs_pp for performing RIXS calculation with $\pi$-$\pi$ polarization and rixs_ps for performing RIXS calculation with $\pi$-$\sigma$ polarization. 2. 2.

Set the control parameters for ed.x, xas.x and rixs.x in file config.in,

File: config.in

&control num_val_orbs=12 ! Number of total valence orbitals (t2g): 2x6=12 num_core_orbs=12 ! Number of total core orbitals (2p): 2x6=12 ed_solver=0 ! Type of ED solver, full diagonalization of H_{i} neval=220 ! Number of eigenvalues obtained nvector=2 ! Number of eigenvectors obtained idump=.true. ! Dump eigenvectors to file eigvec.xxx num_gs=2 ! Numer of gound states used for XAS and RIXS ! calculations linsys_tol=1E-10 ! Tolerance for the termination of solving ! linear equations nkryl=500 ! Maximum number of Krylov vectors when building ! Krylov space gamma_in=2.5 ! Core-hole life-time in eV omega_in=-540.4 ! Incident x-ray energy at which RIXS ! calculations are performed &end 3. 3.

Get parameters of hopping and Coulomb interaction $t_{\alpha,\beta}$ and $U_{\alpha,\beta,\gamma,\delta}$, and write them in files: hopping_i.in, hopping_n.in, coulomb_i.in and coulomb_n.in. 4. 4.

Get Fock basis and write to files: fock_i.in, fock_n.in and fock_f.in. Please note that when setting up the Fock basis, we only consider the valence orbitals because the core orbitals will be explicitly considered within ed.x, xas.x and rixs.x. Here, for $\hat{H}_{i}$, there will be 9 electrons occupying 12 valence orbitals, and for $\hat{H}_{n}$ with a core-hole, there will be 10 electrons occupying 12 valence orbitals. 5. 5.

Get transition operators for XAS and RIXS and write to files: transop_xas.in, transop_rixs_i.in and transop_rixs_f.in. For XAS, we use isotropic polarization. For RIXS, we do calculations for $\pi$-$\pi$ and $\pi$-$\sigma$ polarization. The scattering plane is chosen as $ac$-plane. The incident and scattered angles are $\theta_{\text{in}}=\pi/6$ and $\theta_{\text{out}}=\pi/3$, respectively.

Then, we go to the ed directory and run

$ mpirun -np 4 ed.x > log.dat

to perform the ED calculation. After ED is finished, we copy eigenvectors to xas, rixs_pp and rixs_ps directories,

$ cp eigvec.* ../xas
$ cp eigvec.* ../rixs_pp
$ cp eigvec.* ../rixs_ps

We can now go to xas directory to do XAS calculation and just type

$ mpirun -np 4 xas.x > log.dat

After the XAS calculation is finished, we can get the XAS spectrum by running a Python script in upper directory

$ get_spectrum.py -N 1000 -ommin -560 -ommax -520 -T 300 -G 2.5  \
>  -off 11755.5 -f xas_with_hopp.dat xas/*poles*

where, we get the XAS spectrum near Ir-$L_{3}$ edge among energy range [ommin+off,ommax+off]. Temperature is set to be 300 K and $\Gamma_{c}=2.5$ eV. We can now plot the XAS spectrum, which is shown in Fig. 6(a).

Finally, we go to the rixs_pp and rixs_ps directories to do RIXS calculations. We need manually change omega_in in file config.in to set the incident x-ray energy at which the RIXS are performed. Here, we set it to be -540.4 eV, which is corresponding to the resonant energy 11.215 keV of the Ir-$L_{3}$ edge. We run the command,

$ mpirun -np 4 rixs.x > log.dat

to launch the RIXS calculations. After the calculations are finished, we run the Python script

$ get_spectrum.py -N 1000 -ommin -0.2 -ommax 1.5 -T 300 -G 0.04   \
>  -f rixs_with_hopp.dat rixs_pp/*poles* rixs_ps/*poles*

to get the RIXS spectrum among energy range [-0.2, 1.5]. $\Gamma$ is set to be 0.04 eV. We repeat the above calculations for the case without hoppings between the two Ir sites by setting $t_{1}=t_{2}=0.0$ in file get_inputs.py. The RIXS spectrum is plotted in Fig. 6(b). Without hoppings, we can find three peaks B, C, and F corresponding to the local $d$-$d$ excitations which are mainly determined by Hund’s coupling $J_{H}$ and SOC $\lambda$. After turning on the hoppings, the peak positions of these $d$-$d$ excitations will be shifted and new peaks A, D and E corresponding to dimer excitations appear in the RIXS spectrum. The onset of dimer excitations in RIXS spectrum have been observed in Ba5AlIr2O11 experimentally [56]. Here, we demonstrate that EDRIXS can be helpful to simulate such dimer excitations in materials.

5.3 Anderson impurity model: charge transfer excitations

In this section, we show an example to simulate RIXS spectrum based on Anderson impurity model in an Os compound [57] with nominal $d^{3}$ configuration to get charge transfer excitations. We use a $t_{2g}$ Anderson impurity model with three bath sites per orbitals. The bath sites are obtained from a converged DFT+DMFT calculations. We first go to the examples directory and run the Python script,

$ cd ${EDRIXS_DIR}/examples/cpc/anderson_impurity
$ ./get_inputs.py

to build necessary directories, prepare input files just like what we do in Section 5.2. For Anderson impurity model, we first need to perform ED calculations in subspaces with fixed total occupancy number to determine the total occupancy number of the ground states.

$ cd search_gs
$ echo "mpirun -np 4 ${EDRIXS_DIR}/bin/ed.x" > mpi_cmd.dat
$ ${EDRIXS_DIR}/bin/search_gs_by_N.py -ntot 24 -nimp 6  \
> -N1 2 -N2 22 -mpi_cmd mpi_cmd.dat
$ cd ..

where, the total number of orbitals is 24 and the first 6 are the impurity orbitals. We perform the ED calculations for the total occupancy number ranging from 2 to 22. For these ED calculations, we only need to get rough results, so we use the ED solver with standard Lanczos algorithm without re-orthogonalization by setting ed_solver=1 in file config.in.

After then ED calculations are done, we find that the total occupancy of the global ground state is 15 in the first line of the file results.dat. Then, we run

$ ./get_fock_basis.py -norb 24 -noccu 15

to set up Fock basis. After that, we can go to the ed directory to perform an ED calculation with total occupancy 15. At this time, we use the ED solver with parallel ARPACK library to get more accurate ground states by setting ed_solver=2.

$ cd ed
$ mpirun -np 4 ${EDRIXS_DIR}/bin/ed.x > log.dat
$ cp eigvec.* ../xas
$ cp eigvec.* ../rixs_pp
$ cp eigvec.* ../rixs_ps
$ cd ..

We go to the xas directory to perform a XAS calculation and then get the XAS spectrum by running,

$ get_spectrum.py  -N 1000 -ommin -20 -ommax 10 -T 50 -G 2.5   \
>  -off 10877.2 -f xas.dat xas/*poles*

Finally, we can perform RIXS calculations for polarization $\pi$-$\pi$ and $\pi$-$\sigma$. We first set omega_in to be -6.2, which is corresponding to the resonant energy 10871 eV at Os-$L_{3}$ edge.

$ cd rixs_pp
$ mpirun -np 4 ${EDRIXS_DIR}/bin/rixs.x > log.dat
$ cd ../rixs_ps
$ mpirun -np 4 ${EDRIXS_DIR}/bin/rixs.x > log.dat
$ cd ..

After that, we run

$ ${EDRIXS_DIR}/bin/get_spectrum.py -N 1000 -ommin -0.2 -ommax 7   \
>  -T 50 -G 0.075 -f rixs.dat rixs_pp/*poles* rixs_ps/*poles*

to get the simulated RIXS spectrum.

The simulated XAS and RIXS spectrum are plotted in Fig. 7. In Fig. 7(b), we can see four peaks below 2 eV, which are corresponding to $d$-$d$ excitations in the $t_{2g}$ subspace, and a broad charge transfer peak in the range 5 $\sim$ 7 eV. More details and results of the RIXS simulations based on DFT+DMFT calculations for this compound will be published in another separate paper [58].

6 Summary and future developments

In this paper, we introduce the open source toolkit EDRIXS to simulate RIXS spectra. We explain the basic theory and algorithms for RIXS simulations and the implementation details of EDRIXS code. We also show three examples to demonstrate its usage. EDRIXS is designed as a platform for theoretical simulations of x-ray scattering spectroscopy and it will be a very helpful toolkit for the x-ray scattering community.

The development of EDRIXS code is still in progress. The plans of future development are likely to be along the following directions. More powerful ED solvers based on algorithms such as quantum chemistry (configuration interaction) [59, 60, 61, 62], NRG with non-abelian symmetries [63] will be implemented to diagonalize even larger size of Anderson impurity model. The Python API will be enhanced to provide more powerful and friendly functionalities. A database for Slater integrals of atoms will be provided. Graphical user interface (GUI) can also be implemented for easy use, especially for the experimentalists at beam lines.

7 Acknowledgments

We would like to thank very helpful discussions with Hu Miao and Yongxin Yao. Y.L.W., M.P.M.D. and G.K. were supported by the US Department of energy, Office of Science, Basic Energy Sciences as a part of the Computational Materials Science Program through the Center for Computational Design of Functional Strongly Correlated Materials and Theoretical Spectroscopy. G.F. was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Early Career Award Program under Award No. 1047478. Work at Brookhaven National Laboratory was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-SC0012704.