X-Ray Absorption Spectra Using BSE

by Christian Vorwerk, Caterina Cocchi, & Keith Gilmore for exciting oxygen

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.

2. The N-K edge in cubic boron nitride

Here, we consider the N K-edge in c-BN. As a first step, create a new directory, which you can call, e.g., BN-XANES and move inside it (the symbol $ indicates the shell prompt).

$ mkdir BN-XANES
$ cd BN-XANES

We start with a ground-state calculation of c-BN. Below is an example of the corresponding input file.

<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>

Copy and paste the text of the example above into the file input.xml. Replace the $EXCITINGROOT variable with the correct path by using the command:

$ SETUP-excitingroot.sh

Run the ground-state calculation with the usual command:

$ time exciting_smp

After the ground-state calculation is completed, we inspect the output file 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.38038902

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. Paste the following block 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" 
      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 exciting_smp

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 look at the results contained in the EPSILON folder and execute the following commands:

$ cp EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT singlet-TDA
$ PLOT-files.py -f singlet-TDA  -lx 'Energy [eV]'  -ly 'Im $\varepsilon_M$'  -t 'Macroscopic dielectric function'  -g  -rc  -cy 3  -nl

The result is stored in PLOT.png and PLOT.eps and should resemble this plot:

cBN-XAS-new.png
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 the next example, we will calculate the Ti L2,3-edge XANES of the rutile phase of TiO2. 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 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 core 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 input file for performing the ground-state calculation on rutile TiO2. Copy it into the file input.xml and change the path in speciespath with SETUP-excitingroot.sh.

<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>

Start the calculation by typing the usual command: 
$ time exciting_smp

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.2875257    
 n =  2, l =  0, k =  1 :   -19.33482231    
 n =  2, l =  1, k =  1 :   -16.08151430    
 n =  2, l =  1, k =  2 :   -15.86989435

The core levels we are interested in are 2p states, which have principal quantum number n=2, and angular quantum number l=1. $\kappa$ is 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" 
      vkloff="0.05 0.15 0.25"
      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 exciting_smp

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. Type:

$ cp EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT singlet-TDA
$ PLOT-files.py -f singlet-TDA  -lx 'Energy [eV]'  -ly 'Im $\varepsilon_M$'  -t 'Macroscopic dielectric function'  -x 420 470  -g  -rc  -cy 3  -nl
TiO2-XAS-new.png
Exercise 2
  • 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: Rename the EPSILON directory containing the "L23" results to EPSILON-L23. Repeat the same procedure after the "L2" calculation to produce EPSILON-L2 and after the "L3" calculation to produce EPSILON-L3). How do the three spectra differ? Why?
☛ Solution: L2 and L3 spectra
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License