X-Ray Absorption Spectra Using BSE

for lithium version of Exciting

---

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 the 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. The 1s state corresponds to the K edge in spectroscopic notation.

---

---

0. Preparation

Proceed in the same way as for Excited States from TDDFT.

---

1. Introduction

In order to solve the Bethe-Salpeter equation (BSE), we first need eigenenergies and wave functions from a groundstate calculation as input. In the present methodology the core-edge of interest must be in the valence or otherwise set as a "local orbital", and cannot be a core state. If the corresponding 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 for the N-K edge in cubic BN further down below.

A description of basic setup and parameters can be found in Excited States From TDDFT. Further, all the input parameters are described in Input Reference. Observe that for production results the X-ray absorption spectra need to be converged regarding several parameters such as the size of the k-mesh, see Excited States From TDDFT.

Notice that in any input file in the following the $EXCITINGROOT variable has to be replaced by its actual value (the path itself).

---

2. The Li-K edge in lithium fluoride

Read through the text below in order to complete "Exercise 1"

Let us start with the Li-K edge in LiF. The Li 1s state is quite shallow in energy and it is not put 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"/>
...
    <lorb l="0">
      <wf matchingOrder="0" trialEnergy="0.1500" searchE="false"/>
      <wf matchingOrder="1" trialEnergy="0.1500" searchE="false"/>
      <wf matchingOrder="0" trialEnergy="-1.8034" searchE="true"/>
    </lorb>
...

Note that core="false". This means that we don't need to include any new local orbital, and we can directly use the Li species file $EXCITINGROOT/species/Li.xml.

Below is an example input.xml file which can be used for the LiF groundstate and the Li-K edge calculations xs in the exercise (do not forget to change the path in speciespath).

<input>
 
  <title>Lithium fluoride</title>
 
  <structure speciespath="$EXCITINGROOT/species">
    <crystal>
      <basevect>3.80402 3.80402 0.00000</basevect>
      <basevect>3.80402 0.00000 3.80402</basevect>
      <basevect>0.00000 3.80402 3.80402</basevect>
    </crystal>
    <species speciesfile="Li.xml">
      <atom coord="0.0  0.0  0.0" bfcmt="0.0 0.0 0.0" />
    </species>
    <species speciesfile="F.xml">
      <atom coord="0.5  0.5  0.5" bfcmt="0.0 0.0 0.0" />
    </species>
  </structure>
 
  <groundstate lradstep="2" 
               lmaxvr="8" 
               ngridk="4  4  4" 
               vkloff="0.05 0.15 0.25" />
 
  <xs xstype="BSE" nosym="true" 
      ngridq="2 2 2" ngridk="2 2 2"
      vkloff="0.05 0.15 0.25" nempty="20" 
      lmaxapwwf="3" lmaxemat="3" 
      gqmax="2.0" broad="0.0073499"
      tevout="true" scissor="0.00000">
    <energywindow intv="0.0 5.0" points="2000" />
    <screening screentype="full" nempty="115" />
    <BSE bsetype="singlet" nstlbsemat="1 1 1 15" nstlbse="1 1 1 15" />
    <qpointset>
      <qpoint>0.0 0.0 0.0</qpoint>
    </qpointset>
  </xs>
 
</input>

The attributes nstlbsemat and nstlbse in BSE give the band ranges for the computation of the matrix elements and for diagonalizing the Hamiltionian matrix, respectively. Here, the band ranges are set to "1 1 1 15", where the first part, 1 1, refers to the core-edge state as shown, e.g., in the ground-state output EIGVAL.OUT and the second part, 1 15, gives the conduction states counted from the Fermi level.

Setting bsetype="singlet" in BSE corresponds to including all terms in the effective Hamiltonian (diagonal, direct and exchange). Instead, bsetype="RPA" gives the random-phase approximation (RPA), which does not include the direct term and, therefore, does not give rise to bound core-exciton states.

Important and useful settings for the input file

In order to not recalculate the groundstate

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

To avoid recalculating the screening and matrix elements, set within xs

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

This is useful for instance when computing 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". Note that the present way of performing this setup might change in future versions of the code.

XANES spectra and core-exciton output

The file INFOXS.OUT shows the general messages for the calculation procedure. The absorption spectra is given in EPSILON_BSEsinglet_SCRfull_OC11.OUT, with a Lorentzian broadening broad defined in energywindow inside the xs element. 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 can, e.g., 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

The absorption spectra correspond to the imaginary part of the dielectric function, which can be plotted, e.g., by invoking xmgrace.

$ xmgrace dielectricfunction_Im.agr &

Exercise 1

• Set up and execute a groundstate calculation for LiF (or use result from previous tutorials).
• Perform BSE calculations for the Li-K edge bsetype="singlet".
• Plot the absorption spectra, e.g., using EPSILON_BSEsinglet_SCRfull_OC11.OUT.xml.
• 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".

The figure below shows the Li-K edge spectra in LiF using 15 conduction states for different choice of mesh and gqmax.

---

3. The N-K edge in cubic boron nitride

"Exercise 2" and "Exercise 3" are given at the end of the section

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" bfcmt="0.0 0.0 0.0" />
                </species>
                <species speciesfile="N.xml">
                        <atom coord="0.25 0.25 0.25" bfcmt="0.0 0.0 0.0" />
                </species>
        </structure>
 
        <groundstate gmaxvr="14.0" ngridk="4 4 4"
            swidth="0.01d0" findlinentype="simple" mixer="lin" />
 
</input>

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 include the 1s state of N as a "local orbital". The procedure is shown in the following.

Including a local orbital

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. After successful convergence, check the energy of the state in EVALCORE.OUT, it may look like

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

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

Now, we will make a separate groundstate calculation, including the local orbital. Copy the input.xml file and the $EXCITINGROOT/species/N.xml to a new working directory. It is useful to change the name of the species file to, e.g., N_XANES.xml to not confuse it with others. In 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".

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

In N_XANES.xml, insert also a local orbital by adding

...
<!-- Here below a local orbital is added for 1s -->
  <lorb 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"/>
  </lorb>
...

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

In order to complete the input file, do not forget to change the path in speciespath to the one corresponding to the file N_XANES.xml in the new input.xml. Then, you can run excitingser in the usual way. If calculations with a new local orbital do not converge, it might be necessary to set a different trialEnergy. Also, it might be necessary to explicitly set the parameters deband (increase default value) and epsband (decrease, if needed). See Input Reference for details.

After a successful groundstate calculation including the new local orbital, it is time for considering the XANES computation. An example is given below for bsetype="singlet".

    ...
    <xs xstype="BSE" nosym="true" ngridk="2 2 2" ngridq="2 2 2"
        vkloff="0.05 0.15 0.25"
        nempty="60" lmaxapwwf="5" lmaxemat="5" gqmax="2.0" broad="0.0183747"
        tevout="true">
        <energywindow intv="13.4 16.5" points="1550" />
               <screening screentype="full" nempty="115" />
               <BSE bsetype="singlet" nstlbsemat="1 1 1 15" nstlbse="1 1 1 15" />
        <qpointset>
            <qpoint>0.0 0.0 0.0</qpoint>
        </qpointset>
    </xs>
    ...

Exercise 2

• Perform a regular ground-state calculation for c-BN.
• Include the N 1s local orbital in a new groundstate calculation, compare EIGVAL.OUT.
• Perform BSE calculation for the N-K edge bsetype="singlet".
• Plot the absorption spectra, e.g. using EPSILON_BSEsinglet_SCRfull_OC11.OUT.xml.
• 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 (ngridq="4 4 4" and ngridk="4 4 4".

The figure below shows the B-K edge in c-BN using 15 conduction states and for different choice of mesh and gqmax.

Exercise 3

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

Some tips and hints for production calculations

[show tips]

• Copy atom species files to a new file for including local orbitals.
• If convergence problems appear, check trialEnergy in species file, deband, and epsband in input.xml.
• 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.
• Typical values: lmaxapwwf="5", lmaxemat="5". For gqmax check convergence.
• 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).

---