The RPA and RPA+BSE Dielectric functions

This tutorial will go through how the macroscopic dielectric function ϵM(ω)\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 kk-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 ωmin\omega_{min}; broadening at ωmax\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 ϵM(ω)\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 ϵM\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. LiF absorption spectrum

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


Edit This Page