X-ray Absorption Spectra Using BSE

by Weine Olovsson, Stephan Sagmeister, & Lorenzo Pardini for exciting boron

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 for the Li K-edge in lithium fluoride and N K-edge in cubic boron nitride. We remind that, in spectroscopic notation, the K-edge corresponds to transitions from 1s core state to conduction bands.

---

0. Define relevant environment variables

Before starting, be sure that relevant environment variables are already defined as specified in How to Set Environment Variables for Tutorial Scripts. In particular, the most relevant script for this tutorial is

• PLOT-loss-function.py, a Python script for generating plots of the loss function.

From now on, the symbol $ will indicate the shell prompt.

Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.

---

1. Introduction

In order to solve the Bethe-Salpeter equation (BSE), we first need eigenenergies and wave functions from a ground-state calculation as input. In the present methodology, the core-edge of interest must be included in the valence state or, otherwise, must be set as a "local orbital" and in any case cannot be treated as a core-state. If the relevant core-edge state is not in the valence and not already set as a local orbital, it is possible to do so by hand. This procedure will be shown explicitly in the example further down below for the N-K edge in cubic BN.

A description of basic setup and parameters can be found in Excited states from BSE. Further, 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 (see Excited states from BSE).

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

---

2. The Li-K edge in lithium fluoride

i) Preparation of the input file

Before starting any calculation, create a directory named LiF-XANES and move into it.

$ mkdir LiF-XANES
$ cd LiF-XANES

Let us now start with the Li K-edge in LiF. In this standard notation, the considered core edge is the Li 1s state. This state is quite shallow in energy and exciting does not include it into the core by default. This can be directly seen by checking the Li atom species file $EXCITINGROOT/species/Li.xml

...
      <atomicState n="1" l="0" kappa="1" occ="2.00000" core="false"/>
...

Note that in the element atomicState, the attribute core is set to "false". This means that we do not need to include a new local orbital to describe this state as a valence state, and we can directly use the Li species file $EXCITINGROOT/species/Li.xml.

The file below is an example for the input.xml file which can be used for the LiF groundstate calculation and for the calculations of the the Li-K edge X-ray absorption spectrum (do not forget to change the path in speciespath).

<input>
 
   <title>Lithium fluoride</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="7.608">
         <basevect> 0.5  0.5  0.0 </basevect>
         <basevect> 0.5  0.0  0.5 </basevect>
         <basevect> 0.0  0.5  0.5 </basevect>
      </crystal>
      <species speciesfile="Li.xml">
         <atom coord="0.0  0.0  0.0" />
      </species>
      <species speciesfile="F.xml">
         <atom coord="0.5  0.5  0.5" />
      </species>
   </structure>
 
   <groundstate 
      ngridk="4  4  4" 
      gmaxvr="14" >
   </groundstate>
 
   <xs 
      xstype="BSE" 
      ngridk="2 2 2" 
      vkloff="0.097 0.273 0.493"
      ngridq="2 2 2"
      nempty="30"  
      gqmax="3.0" 
      broad="0.007"
      tevout="true" >
 
      <energywindow 
         intv="1.0 3.0" 
         points="2000" />
 
      <screening 
         screentype="full" 
         nempty="100" />
 
      <BSE 
         bsetype="singlet" 
         nstlbsemat="1 1 1 4" 
         nstlbse="1 1 1 4" />
 
      <qpointset>
         <qpoint>0.0 0.0 0.0</qpoint>
      </qpointset>
 
   </xs>
 
</input>

Elements and attributes appearing in this input file are discussed in Excited states from BSE. In particular, the first two integers of the attributes nstlbse and nstlbsemat are both set to 1, since, in this example only X-ray absorption from the 1s core state is considered.

ii) Calculation of XANES spectra and core-exciton output

To run the calculation, copy and paste the text of the example above into the file input.xml inside the working directory LiF-XANES and start the calculation by invoking the exciting executable:

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

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

iii) Useful settings for the input file

In practical calculations, several runs have to be performed which use the same ground-state results and different settings for the excited-state BSE part. In order to reduce the computational effort, one can avoid to recalculate the ground state. This can be done by adding the do attribute to the groundstate element and set it as follows:

...
   <groundstate 
      do="skip" 
      ... >
...

Furthermore, one can avoid to recalculate screening and matrix elements, which are set inside the element xs. The procedure which allows to skip these consists in adding the subelement plan inside the element xs as follows

...
   <xs>
      ...
      <plan> 
         <doonly task="bse" />
      </plan>
      ...
   </xs>
...

In the example reported above, the execution of the part of the code related to the element screening will be omitted. This is useful when, once the calculation related to the input file reported above has been performed, one wishes to compute a different bsetype afterwards, such as bsetype = "RPA", or when changing the band ranges for the Hamiltonian diagonalization, e.g., nstlbse = "1 1 1 25".

Exercise 1

• Check after bound core-exciton states in EXCITON_BSEsinglet_SCRfull_OC11.OUT.
• Perform a RPA calculation (bsetype = "RPA").
• Plot RPA absorption spectra, e.g., using EPSILON_BSERPA_SCRfull_OC11.OUT.xml and compare with BSE.
• Use a denser mesh for the BSE (singlet) calculations by setting ngridq = "4 4 4" and ngridk = "4 4 4".
• Increase the range of empy states included in the BSE calculation by setting nstlbse = "1 1 1 15" and nstlbsemat = "1 1 1 15"

---

3. The N-K edge in cubic boron nitride

i) Preparation of the input file for the groundstate calculation

Here, we consider the N K-edge in c-BN. Below is an example which can be used for the regular c-BN ground-state calculation (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"
      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.

ii) Including a local orbital

The N 1s state is relatively deep, around 1 keV, and it is by default a core-state in the species file $EXCITINGROOT/species/N.xml. If you look inside this file, you will find the following line.

...
     <atomicState n="1" l="0" kappa="1" occ="2.00000" core="true"/>
...

In order to perform the BSE calculations, we must treat the 1s state of N as a "local orbital". The procedure is shown in the following.

First, we need to estimate the energy location of the state in order to make a good guess for the local orbital trial energy. This is done by a regular groundstate calculation using the default species file $EXCITINGROOT/species/N.xml:

$ time excitingser

After successful convergence, check the energy of the state in EVALCORE.OUT, that should look like

Species :    1 (B), atom :    1
 n =  1, l =  0, k =  1 :   -5.952183887    

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

Now, we will make a separate ground-state calculation, including the local orbital. In order to include the local orbital among the valence states, the species file which describes the atomic configuration of the nitrogen atom must be changed. Create a new working directory, e.g., BN-XANES-lo, and copy the previous input.xml file as well as the files $EXCITINGROOT/species/N.xml and $EXCITINGROOT/species/B.xml in it. It is useful to change the name of the species file to, e.g., N_XANES.xml to avoid confusion among the species files. Inside the new species file N_XANES.xml, remove the core-edge state from the core by changing core = "true" to "false" in the line corresponding to the atomicState n = "1" and l = "0".

...
     <atomicState n="1" l="0" kappa="1" occ="2.00000" core="false"/>
...

Furthermore, in N_XANES.xml, insert also a local orbital in the element basis, which now should look like

...
   <basis>
      <default type="lapw" trialEnergy="0.1500" searchE="false"/>
      <custom l="0" type="apw+lo" trialEnergy="0.1500" searchE="true"/>
      <custom l="1" type="apw+lo" trialEnergy="0.1500" searchE="true"/>
      <!-- Here below a local orbital is added for the state 1s -->
      <lo l="0">
         <wf matchingOrder="0" trialEnergy="  0.1500" searchE="false"/>
         <wf matchingOrder="1" trialEnergy="  0.1500" searchE="false"/>
         <wf matchingOrder="0" trialEnergy="-13.5000" searchE="true"/>
      </lo>
   </basis>        
...

The most important issue is to set a "good enough" trialEnergy in the last row. This value should be a bit lower than the result for the core state. Here, a value of -13.5 Ha is selected (to be compared with -13.38044431 Ha from above).

In order to complete the input file input.xml, two further modifications must be applied:

• You must change the path in speciespath to the one corresponding to the file N_XANES.xml, i.e., to the current directory "./".

...
   <structure speciespath="./">
...

• You must change the name of the species file for the Nytrogen atom to N_XANES.xml

...
      <species speciesfile="N_XANES.xml">
...

Then, you can run excitingser in the usual way. If the calculation with a new local orbital does not converge, it might be necessary to set a different trialEnergy. Also, it might be necessary to explicitly set inside groundstate the attributes deband (by increasing the default value) and epsband (decreasing the default value, if needed). See Input Reference for details.

After a successful ground-state calculation including the new local orbital, it is now time to perform the XANES calculation. To this purpose, add do = "skip" to the element groundstate and then append the following block of instructions after the groundstate:

...
   <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 
         bsetype="singlet" 
         nstlbsemat="1 1 1 15" 
         nstlbse="1 1 1 15" />
 
      <qpointset>
      <qpoint>0.0 0.0 0.0</qpoint>
      </qpointset>
 
   </xs>
...

Then, run the calculation with the usual command

$ time excitingser

Once the calculation is finished, the XANES spectra can be plotted using the procedure explained above in the case of LiF, i.e., running the following commands:

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

Exercise 2

• Check after bound core-exciton states in EXCITON_BSEsinglet_SCRfull_OC11.OUT.
• Perform an RPA calculation (bsetype = "RPA").
• Plot "RPA" absorption spectra, e.g., using EPSILON_BSERPA_SCRfull_OC11.OUT.xml and compare with "BSE".
• Use a denser mesh for the "BSE" (singlet) calculations (ngridq = "4 4 4" and ngridk = "4 4 4").

Exercise 3

• Check the sensitivity of the LiF and c-BN spectra to different parameters, e.g., the number of conduction states (specified by the last integer in nstlbsemat and nstlbse).
• Calculate the boron-K edge in c-BN.

iii) Tips and hints for production calculations

[show tips]

• Investigate the sensitivity of the spectra on the number of conduction states.
• For studying the very near-edge spectra, fewer conduction states are typically needed.
• For comparison with a low resolution experiment, a less dense mesh in xs could be sufficient.
• Check convergence for gqmax.
• In general, BSE computations may be impractical for unit cells containing many atoms.
• Core-exciton binding energies are more difficult to converge than the general spectral shape.
• L_2,3-edge (2p) spectra are more costly than K-edge (more core states included, larger matrix to diagonalize).