q-Dependent BSE Calculations

by Christian Vorwerk & Caterina Cocchi for exciting nitrogen

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

The resulting input file is the following:

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

N. B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT.

$ SETUP-excitingroot.sh

Start now the ground-state SCF calculation and check if it finishes gracefully. 
$ time excitingser &

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. Theoretical background: Electron and x-ray scattering spectroscopy

The inelastic, non-resonant scattering of electrons or x-ray photons with the electrons of a systems are measured in angle-resolved electron energy-loss spectroscopy (AR-EELS) and inelastic X-ray scattering spectroscopy (IXS), respectively. In both spectroscopies, the double differential cross section (DDCS) $\frac{\delta^2 \sigma}{\delta \Omega \omega'}$, which describes the scattering probability of a particle with final energy $\omega'$ into the solid angle $\Omega$. For the scattering of x-ray photons, the DDCS is proportional to the dynamical structure factor $S(\mathbf{q},\Delta \omega)$

(1)
\begin{align} \left.\frac{\delta^2 \sigma}{\delta \Omega \omega'}\right|_{\textrm{x-ray}} \,\propto\; S(\mathbf{q},\Delta \omega)\;, \end{align}

where $\mathbf{q}=\mathbf{k}_i-\mathbf{k}_f$ is the momentum loss from the initial momentum $\mathbf{k}_i$ to the final one $\mathbf{k}_i$, and $\Delta \omega=\omega_i-\omega_f$ is the corresponding energy loss.

For the scattering of electrons, the DDCS is proportional to the electron energy loss function $L(\mathbf{q},\Delta \omega)$:

(2)
\begin{align} \left.\frac{\delta^2 \sigma}{\delta \Omega \omega'}\right|_{\textrm{electron}} \,\propto\; L(\mathbf{q},\Delta \omega)\;, \end{align}

In linear response theory, both the dynamical structure factor and the loss function can be expressed in terms of the macroscopic dielectric function $\varepsilon_{M}(\mathbf{q},\Delta \epsilon)$. We obtain

(3)
\begin{align} S(\mathbf{q},\Delta \omega)=-\frac{1}{\pi}v^{-1}\!(\mathbf{q})\: \Im\!\left[ \frac{1}{\varepsilon(\mathbf{q},\Delta \omega)} \right], \end{align}

where $v(\mathbf{q})$ is the Coulomb potential, and

(4)
\begin{align} L(\mathbf{q},\Delta \omega)=-\Im\!\left[ \frac{1}{\varepsilon(\mathbf{q},\Delta \omega)} \right]. \end{align}

In the theory part of the tutorial Excited States from BSE, the calculation of the macroscopic dielectric function $\varepsilon_M(\omega)$ is discussed in detail. Further in-depth discussion of this theory can be found here.

3. Optical BSE calculation for finite momentum values

i) Preparation of the input file and run

To start the BSE calculations, move to the parent directory LiF-BSE-q and create a new folder which we name BSE. Then, we move into it:

$ cd ..
$ pwd
/home/tutorials/LiF-BSE-q
$ 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 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, also copied from the GS folder.

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.0 0.0 0.0</qpoint>
         <qpoint>0.25 0.0 0.25</qpoint>
         <qpoint>0.5 0.0 0.5</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 is to 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 are now ready to run the BSE calculation:

$ time excitingser &
ii) Analysis of the results

Once the calculation is successfully completed (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, we move to the folder LOSS and execute the following commands:

$ cd LOSS
$ PLOT-loss-function.py --qmt LOSS_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT.xml LOSS_BSE-singlet-TDA-BAR_SCR-full_QMT002.OUT.xml LOSS_BSE-singlet-TDA-BAR_SCR-full_QMT003.OUT.xml
With these commands, the file PLOT.png is created, which looks like the figure given below.
loss-1.png

We notice that, by increasing transferred momentum, the first peak shifts to higher energies and looses oscillator strength. Likewise, the pronounced maximum at about 25 eV, which dominates the spectrum in the optical limit, is quenched at finite momentum transfer.

4. 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 go back of two parent directoriescreate a new sub-directory F-K-edge and move there:

$ cd ../..
$ pwd
/home/tutorials/LiF-BSE-q
$ mkdir F-K-Edge
$ cd F-K-Edge

As usual, we need to copy files from the GS folder into the F-K-Edge 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"
      broad="0.007"
      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.0 0.0 0.0</qpoint>
         <qpoint>0.25 0.0 0.25</qpoint>
         <qpoint>0.5 0.0 0.5</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.
We are now ready to run the BSE calculation: 
$ time excitingser &

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 for the optical spectra: 
$ cd LOSS
$ PLOT-loss-function.py --qmt LOSS_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT.xml LOSS_BSE-singlet-TDA-BAR_SCR-full_QMT002.OUT.xml LOSS_BSE-singlet-TDA-BAR_SCR-full_QMT003.OUT.xml

With these commands, the file PLOT.png is created, which looks like the figure given below.
loss-2.png

In this case, we notice less changes in the loss function compared to the previous case in the optical region. The most relevant changes in the spectrum upon transferred momentum are the emergence of a pre-preak which is dark at $\mathbf{q}=0$, and the increasing intensity of the first strong peak.

5.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.
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License