# The RPA and RPA+BSE Dielectric functions

This tutorial will go through how the macroscopic dielectric function $\epsilon_M(\omega)$ is calculated. The random-phase approximation (RPA) macroscopic dielectric function can be calculated in **lmf** using the **OPTICS** category in the ctrl file; see this page for a detailed description. Whilst the dielectric function in the form produced by **lmf** is useful for analysis there are two crucial effects missing. The first is the effect of local fields, and the second is excitonic (electron-hole interaction) effects. These effects can be included through use of the Bethe-Salpeter equation for the polarisation. See this page for a background to the theory.

### Command Summary

```
$ mkdir lif
$ cd lif
$ nano init.lif # paste in structural information (see below)
$ blm --gw init.lif
$ cp actrl.lif ctrl.lif
$ lmfa ctrl.lif
$ cp basp0.lif basp.lif
$ nano ctrl.lif # change gmax, nit, nkabc, gcutb & gcutx (see below)
$ lmf ctrl.lif
$ echo -1 | lmfgwd ctrl.lif
$ lmgwsc --code2 lif
$ echo 6 | lmfgwd ctrl.lif
$ nano ctrl.lif # Edit OPTICS category
$ lmf --opt:woptmc,rdqp ctrl.lif --rs=1,0 -vnit=1
$ cp optdatac.lif optdatabse
$ nano GWinput # Include parameters for bethesalpeter (see below)
$ echo 0 | bethesalpeter
```

### Preliminaries

**blm**, **lmfa**, **lmf**, **lmfgwd**, **lmgwsc** are all assumed to be in your path. The source code for all Questaal executables can be found here.

You should ensure that the executable **bethesalpeter** is also in your path.

### Tutorial

Start by making a directory **lif**

```
$ mkdir lif && cd lif
```

Structural information about the unit cell is then contained in the file *init.lif*. Create this file and paste the following

```
HEADER Li F (Lithium fluoride)
LATTICE
# SPCGRP=225
# A=4.028 B=4.028 C=4.028 ALPHA=90 BETA=90 GAMMA=90
% const a=4.028
ALAT={a} UNITS=A
PLAT= 0.5000000 0.5000000 0.0000000
0.5000000 0.0000000 0.5000000
0.0000000 0.5000000 0.5000000
SITE
ATOM=Li X= 0.0000000 0.0000000 0.0000000
ATOM=F X= 0.5000000 0.5000000 0.5000000
```

Now run **blm** to produce *ctrl.lif*.

```
$ blm --gw init.lif
$ cp actrl.lif ctrl.lif
```

The –gw switch adjusts certain parameters so that a GW calculation can be performed. Now run **lmfa** to calculate the free-atom density and copy this to *basp.lif*, as this is the file read by **lmf**:

```
$ lmfa ctrl.lif
$ cp basp0.lif basp.lif
```

Since no local orbitals were found we do not need to re-run **lmfa**. The next step then is to edit *ctrl.lif* and pick a suitable $k$-mesh and add the value of **GMAX** suggested by the output of **lmfa**.

```
nit=30
nkabc=5
gmax=9.3
```

One important thing to note for this particular example are the parameters **gcutb** and **gcutx**, these are cut-offs used in the **GW**. Change these to

```
gcutb=3.7 gcutx=3.1
```

**blm** may choose negative values for these parameters, which are not compatible and must be changed. The values given were 3.6 and 3.0 respectively, however, these particular choices for this particular system may causes a segmentation fault when calculating the RPA polarisation needed to calculate the self-energy.

We are now ready to run our self consistent LDA calculation. Type the following

```
$ lmf ctrl.lif
```

This should produce a band gap of about 9.4 eV, which is small compared with experiment ~14.2 eV.

#### QSGW calculation

We now want to perform a **QSGW** calculation. To build the input files required by **QSGW**, type the following

```
$ echo -1 | lmfgwd ctrl.lif
```

The main input is *GWinput*. This file, as it stands, is sufficient for this particular example. We can now perform a QSGW calculation on this system:

```
$ lmgwsc --code2 lif
```

Convergence should be reached after 5 iterations (iteration number 4, since counting starts at 0). The gap produced now should be about 16.4~16.5 eV. This is quite high, but can – however – be corrected slightly with a finer k-mesh.

Before we run the **bethesalpeter** program we need to use **lmf** to produce the optical matrix elements, see this page for a description of how the optical matrix elements enter in to the expression for the macroscopic dielectric function. To do this we need to ensure the k-mesh used by **bethesalpeter** is consistent with that used by **lmf**. To do this run the following command

```
$ echo 6 | lmfgwd ctrl.lif
```

We then need to include information in *ctrl.lif* about which pairs of states etc to calculate the optical matrix elements for. We do this through the OPTICS category, see this page

```
OPTICS
MODE=1
LTET=0
WINDOW=0.0 2.0
NPTS=1001
FILBND=1 4
EMPBND=5 50
MEFAC=0
```

MODE, LTET, WINDOW and NPTS do not affect the calculation of the matrix elements. FILBND and EMPBND include the upper and lower bounds on the valence and conduction bands for which we calculate the matrix elements and MEFAC determines how the contribution to these matrix elements from the non-local part of the potential is included. We do not include this correction here (it is handled in **bethesalpeter**, see below) hence MEFAC=0.

Next we run

```
$ lmf --opt:woptmc,rdqp ctrl.lif --rs=1,0 -vnit=1
$ cp optdatac.lif optdatabse
```

where we copy the file produced to that read in by **bethesalpeter**. Before we perform the Bethe-Salpeter calculation we need to include some parameters in *GWinput*.

```
$ nano GWinput
```

The parameters for the calculation are

```
MaxOmega_BSE 2.0
MinOmega_BSE 0.0
numOmega_BSE 8000
NumValStates 4
NumCondStates 4
NLF_corr_BSE on
Qvec_bse 1 1 1 !q vector for BSE
ImagEner1 -0.083 !0.00160735
ImagEner2 0.107 !0.0138375
RestartwithKernel off
RestartwithDiagonal off
```

which are (in order of appearance): the maximum $\omega$ (in Rydberg) for which the dielectric function is calculated; the minumum $\omega$; the number of frequencies; the number of valence states to include in the Bethe-Salpeter equation; the number of conduction states included; LDA energies ARE used to construct the optical matrix elements accounting for a non-local potential; the direction of the perturbing field; Lorentzian broadening (Rydberg) at $\omega_{min}$; broadening at $\omega_{max}$ (broadening can increase linearly); do not read the kernel from a previous calculation; and do not read the diagonalized 2-particle Hamiltonian from a previous calculation. Note that NumValStates and/or NumCondStates can be less than the number of states included in FILBND and EMPBND in *ctrl.lif*. The valence and conduction states used in the Bethe-Salpeter equation are those closest to the Fermi level. The rest of the states are included in $\epsilon_M(\omega)$, but included at RPA level. The reason for doing this is the computational expense of solving the Bethe-Salpeter equation.

We are now ready to run a Bethe-Salpeter equation calculation to get the macroscopic dielectric function. To do this type

```
$ echo 0 | bethesalpeter > out.bse
$ cp eps_BSE.out eps_BSE.outBSE
```

*eps_BSE.out* contains parameters used, followed by three columns. The energy (for numOmegaBSE energies from MinOmega_BSE to MaxOmega_BSE) in eV are in column one followed by the real and imaginary parts of the macroscopic dielectric function. The imaginary part determines the absorption. We can vary the broadenings to match experiment. Ensure RestartwithKernel and RestartwithDiagonal are switched on to avoid having to calculate and diagonalize the matrix all over again.

Finally, the reason we copied the output file is so that we can compare this with an RPA $\epsilon_M$. To do this type

```
$ echo -1 | bethesalpeter
```

The two outputs can be plotted in a plotting program such as gnuplot. The output from the two calculations above are shown below in the image, along with the experimental spectrum.

Adjust the broadenings to match experiment with RestartwithKernel and RestartwithDiagonal switched on.

Edit This Page