X-ray Absorption Spectra Using BSE

by Christian Vorwerk & Caterina Cocchi for exciting carbon

Purpose: In this tutorial, we present an example of the calculation of the X-ray absorption near-edge structure (XANES) spectra of core-states by solving the Bethe-Salpeter equation (BSE). The general way to set up BSE computations for XANES in exciting is shown, together with examples of the N K-edge in cubic boron-nitride, as well as of the Ti L2,3-edge in rutile titanium dioxide.

1. Introduction

The approach to calculate core spectra from the solution of the Bethe-Salpeter equation (BSE) in exciting is very similar to the approach to optical spectra described in the tutorial Excited states from BSE. The main difference resides in the initial states of the transitions, which are core states. As such, they are solutions of the radial Dirac equation in the muffin-tin spheres.

A description of the basic setup for BSE calculations can be found in Excited states from BSE. In addition, all the input parameters are described in Input Reference. Note that, in order to obtain reliable X-ray absorption spectra, the calculation needs to be converged with respect to several parameters, such as the size of the k-mesh.

Notice that, in any input file in the following, the $EXCITINGROOT variable has to be replaced by its actual value (the path itself) or it can be automatically set by using the command: $ SETUP-excitingroot.sh


From now on, the symbol $indicates the shell prompt. 2. The N-K edge in cubic boron nitride Here, we consider the N K-edge in c-BN. As a first step, we need to compute the ground state of this system, as a starting point for any excited-state calculation. Below is an example of the input file that can be used for the ground-state calculation of c-BN (do not forget to change the path in speciespath). <input> <title>Cubic boron nitride </title> <structure autormt="true" speciespath="$EXCITINGROOT/species">
<crystal scale="6.83136">
<basevect>0.5 0.0 0.5</basevect>
<basevect>0.0 0.5 0.5</basevect>
<basevect>0.5 0.5 0.0</basevect>
</crystal>
<species speciesfile="B.xml" >
<atom coord="0.00 0.00 0.00" />
</species>
<species speciesfile="N.xml">
<atom coord="0.25 0.25 0.25" />
</species>
</structure>

<groundstate
ngridk="4 4 4"
xctype="GGA_PBE_SOL"
gmaxvr="14.0" />

</input>


Go to the parent directory and create a new directory, which you can call, e.g., BN-XANES and move inside it.

$cd ..$ mkdir BN-XANES
$cd BN-XANES  Here, copy and paste the text of the example above for the ground state of BN into the file input.xml. We will now first execute the ground-state calculation with the usual command: $ time excitingser &


After the ground-state calculation is completed, we can inspect one of the output files, namely EVALCORE.OUT, which contains the information about core electrons. Since we are interested here in the excitations from the Nitrogen K-edge, we have to look specifically for the N 1s electron. For this purpose we have to consider the core state with quantum numbers n=1 and l=0:

Species :    2 (N), atom :    1
n =  1, l =  0, k =  1 :   -13.38038416


From the file EVALCORE.OUT we obtain the information that the N 1s electron in cubic BN has an energy of about -13.38 Ha. With this knowledge, we can add to the input file the parameters required to perform a XANES calculation. Since we are dealing with excited-state properties, we have to insert the xs element. The following block can be pasted into the input file:

...
<xs
xstype="BSE"
ngridk="2 2 2"
vkloff="0.05 0.15 0.25"
ngridq="2 2 2"
nempty="30"
gqmax="3.0"
tevout="true" >

<energywindow
intv="13.4 16.5"
points="1500" />

<screening
screentype="full"
nempty="100" />

<BSE
xas="true"
xasspecies="2"
xasatom="1"
xasedge="K"
bsetype="singlet"
nstlxas="1 15" />

<qpointset>
<qpoint>0.0 0.0 0.0</qpoint>
</qpointset>

</xs>
...


You should be already familiar with many parameters shown above. However, it is useful to highlight a few points here:

• The element energywindow is now chosen from 13.4 to 16.5 Ha. This is consistent with the energy of the N 1s state, which represents the investigated absorption edge;
• The attribute xas = "true" in the BSE element triggers the core-level calculation;
• The attributes xasspecies, xasatom, and xasedge uniquely define the initial states of the transition. In our case, we choose xasspecies = "2" to select Nitrogen, xasatom = "1" to identify the one and only N atom in the unit cell, and the xasedge = "K" to specify that we are performing a K-edge calculation.
• The attribute nstlxas describes the number of unoccupied states used in the BSE-Hamiltonian counting from the lowest unoccupied band. Setting it to nstlxas = "1 15", we use the lowest 15 unoccupied bands in the calculation.

We are now ready to run the calculation with the usual command:

$time excitingser &  As output, many files are generated. Among them: • The file INFOXS.OUT shows the general messages for the calculation procedure. • The absorption spectrum is given in EPSILON_BSEsinglet_SCRfull_OC11.OUT, with a Lorentzian broadening broad defined inside the element xs. • The energy and oscillation strength of the excited states are listed in EXCITON_BSEsinglet_SCRfull_OC11.OUT, where negative energies correspond to bound core-excitons. The XANES spectra, i.e., the absorption spectra corresponding to the imaginary part of the dielectric function, can be plotted using an xml visualization template, generating agr files readable by xmgrace. $ cp EPSILON_BSEsinglet_SCRfull_OC11.OUT.xml tmp.xml
$xsltproc --stringparam filename dielectricfunction$EXCITINGVISUAL/xmldielectricfunction2agr.xsl tmp.xml
$xmgrace dielectricfunction_Im.agr >/dev/null 2>&1 &  The resulting spectrum is displayed below. Exercise 1 • Calculate the N K-edge spectra in the random-phase approximation (RPA) by setting the parameter bsetype = "RPA". Compare the spectra of the BSE and the RPA calculation. • Calculate the boron K-edge in c-BN. 3. The Ti L2,3-edge in rutile TiO2 i) Preparation of the input file and ground-state calculation Before starting any calculation, move to the parent directory, create a directory named TiO2-XANES, and move into it. $ cd ../
$mkdir TiO2-XANES$ cd TiO2-XANES


Let us now start with the Ti L2,3-edge in rutile TiO2. In the standard notation, this means that we consider transitions from the Ti 2p states. First, we need to make sure that the initial state is treated as a core state in the calculation. Therefore, we give a look at the species file $EXCITINGROOT/species/Ti.xml ... <atomicState n="2" l="1" kappa="1" occ="2.00000" core="true"/> <atomicState n="2" l="1" kappa="2" occ="4.00000" core="true"/> ...  Note that in both elements atomicState, the attribute core is set to "true". The two different atomic states correspond to the Ti 2p1/2 and Ti 2 p3/2 states, respectively. The file below is an example of the input.xml file for performing the ground-state calculation on rutile TiO2 (do not forget to change the path in speciespath): <input> <title>TiO2_rutile</title> <structure speciespath="$EXCITINGROOT/species">
<crystal>
<basevect>   8.763555397       0.000000000       0.000000000</basevect>
<basevect>   0.000000000       8.763555397       0.000000000</basevect>
<basevect>   0.000000000       0.000000000       5.610015465</basevect>
</crystal>
<species speciesfile="Ti.xml" rmt="1.8000">
<atom coord="0.0000000000      0.0000000000      0.0000000000"/>
<atom coord="0.5000000000      0.5000000000      0.5000000000"/>
</species>
<species speciesfile="O.xml" rmt="1.8000">
<atom coord="0.3050853616      0.3050853616      0.0000000000"/>
<atom coord="0.6949146384      0.6949146384      0.0000000000"/>
<atom coord="0.1949146384      0.8050853616      0.5000000000"/>
<atom coord="0.8050853616      0.1949146384      0.5000000000"/>
</species>
</structure>
<groundstate
rgkmax="7.0"
ngridk="4 4 6"
xctype="GGA_PBE_SOL"
epsocc="1.0d-6">
</groundstate>
</input>


We start the calculation by typing the usual command:
$time excitingser &  Once the calculation is finished, we can identify the Ti 2p states in the output file EVALCORE.OUT, and with them determine the energywindow for the XANES calculation. To do so, we open EVALCORE.OUT and look for the following lines: Species : 1 (Ti), atom : 1 n = 1, l = 0, k = 1 : -178.2875229 n = 2, l = 0, k = 1 : -19.33482209 n = 2, l = 1, k = 1 : -16.08151432 n = 2, l = 1, k = 2 : -15.86989436  The core levels we are interested in are 2p states. Hence, they have principal quantum number n=2, and angular quantum number l=1.$\kappais a unique combination of the quantum numbers l and J given as (1) \begin{align} \kappa=\begin{cases}-l-1 & \text{for } J=l+\frac{1}{2}\\l & \text{for } J=l-\frac{1}{2} \end{cases}. \end{align} The states we are interested in are Ti 2p1/2 (\kappa$=1) at about -16.1 Ha and Ti 2p3/2 ($\kappa$=-2) at about -15.87 Ha. Thus, we find the Ti 2p states at energies between -16.08 and -15.87 Ha. To make sure that our energy window is big enough to also include excitonic effects, we choose energywindow = "15.0 18.0". ii) Calculation of XANES spectra and core-exciton output Now we have to add the xs element to the input file: ... <xs xstype="BSE" ngridk="2 2 3" ngridq="2 2 3" rgkmax="7.0" lmaxapw="12" vkloff="0.05 0.15 0.25" reduceq="false" reducek="false" nempty="30" gqmax="1.0" broad="0.018" tevout="true" > <energywindow intv="15.0 18.0" points="1500" /> <screening screentype="full" nempty="75" /> <BSE xas="true" xasspecies="1" xasatom="2" xasedge="L23" bsetype="singlet" nstlxas="1 20" /> <qpointset> <qpoint>0.0 0.0 0.0</qpoint> </qpointset> </xs> ...  We can also set the attribute do = "skip" in the element groundstate to avoid an unnecessary recalculation of the groundstate. ... <groundstate do="skip" ... > ...  To run the calculation, invoke the exciting executable: $ time excitingser


The XANES spectra, i.e., the absorption spectra corresponding to the imaginary part of the dielectric function, can be plotted using an xml visualization template, generating agr files readable by xmgrace.

$cp EPSILON_BSEsinglet_SCRfull_OC11.OUT.xml tmp.xml$ xsltproc --stringparam filename dielectricfunction $EXCITINGVISUAL/xmldielectricfunction2agr.xsl tmp.xml$ xmgrace dielectricfunction_Im.agr >/dev/null 2>&1 &


Notice, that the energy scale is given in units of eV.

Exercise 1
• Check after bound core-exciton states in EXCITON_BSEsinglet_SCRfull_OC11.OUT.
• Calculate the L2 and L3 subedges by setting the xasedge attribute to "L2" and "L3", respectively. Plot the different spectra in one plot.
• Perform a RPA calculation (bsetype = "RPA").
• Plot RPA absorption spectra, e.g., using EPSILON_BSERPA_SCRfull_OC11.OUT.xml and compare with BSE.