X-Ray Absorption Spectra Using BSE

by Christian Vorwerk & Caterina Cocchi for exciting nitrogen

Purpose: In this tutorial, we present an example of the calculation of the X-ray absorption near-edge structure (XANES) spectra by solving the Bethe-Salpeter equation (BSE). For this purpose, we show the 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>


Create a new directory, which you can call, e.g., BN-XANES and move inside it.

$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.38038417  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 (do not forget to change the path in speciespath): ... <xs xstype="BSE" ngridk="2 2 2" vkloff="0.05 0.15 0.25" ngridq="2 2 2" nempty="30" gqmax="3.0" broad="0.018" 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 (do not forget to skip the ground-state run by inserting the attribute do = "skip" inside the groundstate element): $ time excitingser &


The calculation will run for a few seconds. You can follow the progress on the shell, as described in the Tutorial Excited states from BSE. The produced output is organized as in BSE calculations run in the optical region (see Excited states from BSE). A number of files are printed in the working directory. Among them, the file INFOXS.OUT shows the general messages for the calculation procedure. The physically relevant output quantities are included in the sub-folders named EPSILON, EXCITON, and LOSS. Since we are interested in the X-ray absorption spectrum of cubic-BN, given by the imaginary part of the dielectric function, we move into the EPSILON folder and execute the following commands:

$cd EPSILON$ cp EPSILON_BSE-singlet-TDA-BAR_SCR-full_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 As next example, we will calculate the XANES of the rutile phase of TiO2 from the Ti L2,3-edge. Before starting the new calculation, go back to the parent directory, create a directory named TiO2-XANES, and move into it: $ cd ../../
$mkdir TiO2-XANES$ cd TiO2-XANES


In the standard notation, exciting the Ti L2,3-edge in rutile TiO2 corresponds to considering transitions from all the Ti 2p states. First, we need to make sure that the initial states are treated as core states in exciting. Therefore, we give a look at the species file $EXCITINGROOT/species/Ti.xml. By typing $ cat $EXCITINGROOT/species/Ti.xml  you will see on the screen a list of atomic states. We are interested in the first few lines. The 2p-electrons, given by ... <atomicState n="2" l="1" kappa="1" occ="2.00000" core="true"/> <atomicState n="2" l="1" kappa="2" occ="4.00000" core="true"/> ...  are indeed treated as states. In fact, in both lines 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.08151433 n = 2, l = 1, k = 2 : -15.86989437  The core levels we are interested in are 2p states, which have principal quantum number n=2, and angular quantum number l=1.$\kappais a unique combination of the quantum numbers l and J given by (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 Ti 2p1/2 states (\kappa$=1) are at about -16.1 Ha and the Ti 2p3/2 ones ($\kappa$=-2) are at about -15.87 Ha. Thus, we define the energy window for our XANES by setting energywindow = "15.0 18.0". ##### ii) Calculation of XANES spectra and core-exciton output Now we can start the XANES calculation by adding 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> ...  Remember to set the attribute do = "skip" inside the element groundstate to skip the ground-state calculation. To run the calculation, invoke the exciting executable: $ time excitingser &


Once the calculation is finished (it will take slightly longer than the one on c-BN above), we can inspect the computed XANES spectra. To do so, as usual, we consider the imaginary part of the dielectric function, which is included in the EPSILON folder. By proceeding as before we type:

$cd EPSILON$ cp EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT.xml tmp.xml
$xsltproc --stringparam filename dielectricfunction$EXCITINGVISUAL/xmldielectricfunction2agr.xsl tmp.xml
\$ xmgrace dielectricfunction_Im.agr >/dev/null 2>&1 &

##### Exercise
• Perform a RPA calculation (bsetype = "RPA").
• Plot RPA absorption spectra, e.g., using EPSILON_BSE-RPA-TDA-BAR_SCR-full_OC11.OUT.xml inside the EPSILON subfolder and compare with the spectrum previously calculated from the full BSE. How do the spectra differ? Why?
• Search for bound core-exciton states in the file EXCITON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT inside the EXCITON subfolder.
• Calculate the L2 and L3 subedges by setting the xasedge attribute to "L2" and "L3", respectively. Plot the different spectra in one plot. (Suggestion: perform the calculation for each sub-edge in a separate sub-folder with respect to the working directory TiO2-XANES). How do the three spectra differ? Why?