q-Dependent BSE calculations

by Christian Vorwerk & Caterina Cocchi for exciting oxygen

Purpose: In this tutorial you will learn how to perform Bethe-Salpeter equation (BSE) calculations to obtain x-ray and electron scattering spectra for finite momentum transfer. As an example, the valence and core scattering spectra of LiF are calculated.

### 1. Preliminary step: Ground-state calculation

Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts.

Important note: All input parameters are given in atomic units!

As a preliminary step to calculate excited-state properties from BSE, a ground-state calculation has to be performed. In this tutorial we consider as an example LiF. Create a directory named LiF-BSE-q and move into it.

$mkdir LiF-BSE-q$ cd LiF-BSE-q


Inside the directory LiF-BSE-q create a sub-directory GS where the preliminary ground state calculation is performed:

$mkdir GS$ cd GS


Inside the GS sub-directory we create the input file for LiF. In the structure element we include the lattice parameter and basis vectors of LiF, which has a rock-salt cubic lattice, as well as the positions of the Li and F atoms. In the groundstate element, we include a 10$\times$10$\times$10 k-point mesh (ngridk) and a value of 14.0 for gmaxvr. This value, which is larger than the default, is needed in view of the excited-state calculation.

Create a new file input.xml and copy the following into it:

<input>

<title>LiF-BSE</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.0000 0.0000 0.0000" /> </species> <species speciesfile="F.xml"> <atom coord="0.5000 0.5000 0.5000" /> </species> </structure> <groundstate do="fromscratch" ngridk="10 10 10" xctype="GGA_PBE_SOL" gmaxvr="14.0"/> </input>  Set the speciespath in the usual way $ SETUP-excitingroot.sh


Now start the ground-state SCF calculation and check that it finishes gracefully.
$time exciting_smp  In case of a successful run the files STATE.OUT and EFERMI.OUT should be present in the directory. These two files are needed as a starting point for the BSE calculation. ### 2. Optical BSE calculation for finite momentum values To start the BSE calculations, move to the parent directory LiF-BSE-q and create a new folder with the name BSE and move into it: $ cd ..
$mkdir BSE$ cd BSE


As anticipated, we need the files STATE.OUT and EFERMI.OUT obtained from the ground state calculation to be inside the new directory. Copy these two files, as well as the input file, from the GS folder into the BSE one.

$cp ../GS/{STATE.OUT,EFERMI.OUT,input.xml} ./  In this way, we can skip the ground-state calculation by setting do = "skip" in the groundstate element of the input file. To perform an excited-state calculation we must include the xs element in the input file: ... <xs xstype="BSE" ngridk="4 4 4" vkloff="0.097 0.273 0.493" ngridq="4 4 4" nempty="30" gqmax="3.0" broad="0.007" scissor="0.20947" tevout="true"> <energywindow intv="0.3 1.0" points="1200"/> <screening screentype="full" nempty="100"/> <BSE bsetype="singlet" nstlbse="1 5 1 4" iqmtrange="1 3"/> <qpointset> <qpoint>0.00 0.0 0.00</qpoint> <qpoint>0.25 0.0 0.25</qpoint> <qpoint>0.50 0.0 0.50</qpoint> </qpointset> </xs> ...  The structure of the xs element is nearly identical to the one described in the tutorial Excited States from BSE. The major changes are the following: • The attribute iqmtrange in the BSE element defines the range of momentum transfer vectors$\mathbf{q}$for which the calculation should be performed. The range has to lie within the boundaries of the the element qpointset. • In the qpointset element, more qpoint subelements are added, which define momentum transfer vectors for the calculations. Note that the$\mathbf{q}$-vectors are provided in units of the reciprocal lattice vectors! We can now run the BSE calculation: $ time exciting_smp


Once the calculation completes successfully (it will take a few minutes), a number of files and some folders will be generated in the working directory. Here, we are mostly interested in the calculated loss functions and dynamical structure factors, which can be found in the directory LOSS. The directory contains several files ending with _QMT00x.OUT, where x is the index of the momentum transfer vector. For q=0, the loss function is a 3$\times$3 tensor, and the optical components $L_{xy}(\mathbf{q}=0,\Delta \omega)$ are written to the file ending with _OCxy.OUT.

To plot the scattering spectra, follow the next steps.

• Execute the commands:
$cp LOSS/LOSS_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT LOSS/QMT_1$ cp LOSS/LOSS_BSE-singlet-TDA-BAR_SCR-full_QMT002.OUT LOSS/QMT_2
$cp LOSS/LOSS_BSE-singlet-TDA-BAR_SCR-full_QMT003.OUT LOSS/QMT_3  $ PLOT-files.py -d LOSS  -f QMT_1 QMT_2 QMT_3  -lx 'Energy [eV]'  -ly 'Loss function'  -x 8 27  -y 0 8  -lp 2


These commands create the file PLOT.png:

Notice that increasing the transferred momentum shifts the first peak to higher energies and reduces the oscillator strength. The pronounced maximum above 25 eV, which dominates the spectrum in the optical limit, is quenched at finite momentum transfer.

### 3. X-ray BSE calculation for finite momentum values

Equivalently to the calculation of the optical loss function above, we will now calculate the core loss function for excitations from the F 1s states, i.e., the F K edge of LiF. It is assumed that you are already familiar with this type of excited-state calculations, having performed the tutorial X-ray Absorption Spectra Using BSE. We move back to the parent directory LiF-BSE-q directory and create a new sub-directory NRIXS then move there:

$cd ..$ mkdir NRIXS
$cd NRIXS  As usual, we need to copy files from the GS folder into the NRIXS one, to avoid repeating the ground-state calculation. $ cp ../GS/{STATE.OUT,EFERMI.OUT,input.xml} ./


The ground-state calculation is skipped by setting do = "skip" in the groundstate element of the input file, also copied from the GS folder.
Then, we add the following xs subelement into the input file:
...
<xs
xstype="BSE"
ngridk="4 4 4"
vkloff="0.097 0.273 0.493"
ngridq="4 4 4"
nempty="30"
gqmax="3.0"
scissor="0.20947"
tevout="true">

<energywindow
intv="24.0 24.8"
points="1200"/>

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

<BSE
bsetype="singlet"
xas="true"
xasspecies="2"
xasatom="1"
nstlxas="1 4"
iqmtrange="1 3"/>

<qpointset>
<qpoint>0.00 0.0 0.00</qpoint>
<qpoint>0.25 0.0 0.25</qpoint>
<qpoint>0.50 0.0 0.50</qpoint>
</qpointset>

</xs>
...


The input file is nearly identical to the one used for optical calculations above. We have simply added the attribute xas to trigger a core-level calculation, and we have specified which core states we want to excite using the attributes xasspecies, xasatom, and xasedge. The range of unoccupied states is now defined by the attribute nstlxas. Note that we also have to modify the attribute energywindow, since core excitations occur at much higher energies than the optical ones.

Run the BSE calculation with the usual command:

$time exciting_smp  After the run is successfully completed (it will take a few minutes), the usual output files and sub-folders will appear in the working directory. Again, we are interested in the loss function, which is found inside the folder LOSS. The calculated loss function can be visualized with the same script used above: $ cp LOSS/LOSS_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT   LOSS/QMT_1
$cp LOSS/LOSS_BSE-singlet-TDA-BAR_SCR-full_QMT002.OUT LOSS/QMT_2$ cp LOSS/LOSS_BSE-singlet-TDA-BAR_SCR-full_QMT003.OUT LOSS/QMT_3
$PLOT-files.py -d LOSS -f QMT_1 QMT_2 QMT_3 -lx 'Energy [eV]' -ly 'Loss function' -x 655 675  These commands create the file PLOT.png: We observe less dependence of the loss function on momentum transfer than we saw in the previous case for the optical region. The most relevant changes in the spectrum with increasing momentum transfer are the emergence of a pre-peak that is dark at$\mathbf{q}=0\$ and the increasing intensity of the first strong peak.

### 4. Exercise

• Modify the attribute bsetype to calculate the independent-particle (IP), random-phase approximation RPA and triplet loss function spectrum in the optical range for LiF. How is the spectrum modified? If you have already computed the singlet spectra, you only need to execute the last task of the BSE program flow, i.e., "bse", when switching the bsetype.