# Including Ladder Diagrams in W

This tutorial will go through how to run a QS**GW** calculation whilst including ladder diagrams (vertex contributions) in the screened Coulomb interaction **W**. By default, Questaal will perform a QS**GW** calculation with the polarization calculated at the level of the random phase approximation (RPA). This is calculated in the program hx0fp0 (hx0fp0_sc). To include ladders in **W**, the RPA **W** is first calculated in hx0fp0 and this is then used in the Bethe-Salpeter equation (BSE) for the full polarization matrix (See this page for a background to the theory). The BSE for the improved **W** is calculated by the program bethesalpeter, which is included in the QS**GW** calculation with the switch **UseBSEW** set to 1 in **GWinput** (see below). The new **W** is then used instead of the RPA **W** in hsfp0 (hsfp0_sc) to calculate the screened part of the self-energy.

### Command Summary

```
mkdir NiO
cd NiO
cp path-to-questaal/lm/gwd/test/nio.code/{ctrl.nio,rst.nio,sigm.in,switches-for-lm,basp.nio,site.nio,EPS0001.dat} ./
mv sigm.in sigm
ln -s sigm sigm.nio
lmf -vnit=1 --iactiv=no `cat switches-for-lm` ctrl.nio > out.lmf_QSGW
lmchk ctrl.nio `cat switches-for-lm` --syml~mq~lbl=XGLUT~q=0,1/2,1/2,0,0,0,0,1/2,0,11/32,26/32,11/32,1/2,1/2,1/2
lmf -vnit=1 --rs=1,0 `cat switches-for-lm` --band:fn=syml ctrl.nio
cp bnds.nio bnds.qsgw
lmf -vnit=1 --rs=1,0 -vloptic=1 -vmetal=3 `cat switches-for-lm` ctrl.nio # produces RPA absorption (see below)
echo -1 | lmfgwd `cat switches-for-lm` ctrl.nio
nano GWinput # Change UseBSEW, NumCondStates and add NumValStates (see below)
lmgwsc --iter=1 --maxit=1 --code2 --openmp=N --sym --getsigp --insul=22 nio > out.lmgwsc_bse
lmf -vnit=12 --iactiv=no `cat switches-for-lm` ctrl.nio > out.lmf_QSGWBSE
lmf -vnit=1 --rs=1,0 `cat switches-for-lm` --band:fn=syml ctrl.nio
cp bnds.nio bnds.bse
echo -9,16 / | plbnds -ef=0 -scl=13.6058 -fplot -lbl=X,G,L,U,T bnds.nio
fplot -f plot.plbnds
```

### Preliminaries

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

```
$ mkdir NiO && cd NiO
```

We start this calculation from a converged QS**GW** calculation; some of the steps in this tutorial wont be necessary if you performed the calculation from scratch and these will be highlighted. The screening **W** in the converged calculation was calculated in the RPA, hence the large band gap overestimation (see below). The *ctrl.nio*, *site.nio*, *basp.nio* and *switches-for-lm* files are all located in **~path-to-questaal/lm/gwd/test/nio.code2/**. *switches-for-lm* contains command line arguments that can be changed such as the number of **k**-points, the maximum number of iterations, etc. The QS**GW** self-energy that is read by lmf is also contained in this directory and called *sigm.in*. Finally the file *rst.nio* that is needed for this calculation is contained in this directory. These files would all be produced if one had started this calculation from scratch (See this page on how to run a basic QS**GW** calculation for Si). The sigm file that is read by **lmf** is called sigm.nio and, again, this file would have been created by the script **lmgwsc**.

```
$ cp ~path-to-questaal/lm/gwd/test/nio.code2/{ctrl.nio,rst.nio,site.nio,switches-for-lm,basp.nio,sigm.in,EPS0001.dat} ./
```

The *EPS0001.dat* file contains the macroscopic dielectric function (columns 5 and 6 for the real and imaginary parts) as a function of energy in Rydbergs (column 4). Note that it may be necessary for plotting purposes to run the command

```
$ sed -i -e 's/D/E/g' EPS0001.dat
```

so that exponents in this file are written with an E instead of D. We copy this file to compare with other methods for calculating the macroscopic dielectric function below.

The file *sigm.in* created in **lmgwsc** would actually be called *sigm.nio* and **lmgwsc** would have created a link to *sigm*. To do this here we type

```
$ mv sigm.in sigm
$ ln -s sigm sigm.nio
```

Now run **lmf** to produce the H_{0} and neccessary files needed for **GW**

```
$ lmf --iactiv=no -vnit=12 `cat switches-for-lm` ctrl.nio > out.lmf_QSGW
```

The –iactiv=no turns off the interactive mode and -vnit=12 means the calculation will perform a maximum of 12 SCF iterations. Note that if these switches are placed before the `catswitches-for-lm`

command then these values will take precedence over the values in *switches-for-lm*, however, if the are placed after then the values in the file will be used.

If we grep for the bandgap (grep gap *), we should find that **lmf** produces a bandgap of 4.86 eV, an overestimate of 0.5~0.6eV with respect to experiment. A QS**GW** band structure can be produced following the usual steps, which are discussed below.

We will now create the file needed for plotting the band structure. To produce the *syml.nio* file containing the high-symmetry points and paths around the Brillouin zone, run the command

```
$ lmchk ctrl.nio `cat switches-for-lm` --syml~mq~lbl=XGLUT~q=0,1/2,1/2,0,0,0,0,1/2,0,11/32,26/32,11/32,1/2,1/2,1/2
```

then produce the band structure in the usual way

```
$ lmf -vnit=1 --rs=1,0 `cat switches-for-lm` --band:fn=syml ctrl.nio
$ cp bnds.nio bnds.qsgw
```

We can also produce an RPA optical absorption spectrum here using **lmf**. The necessary categories should be in the *ctrl.nio* file. To do this run

```
$ lmf --iactiv=no -vnit=1 --rs=1,0 -vmetal=3 -vloptic=1 `cat switches-for-lm` ctrl.nio
```

The extra switches are discussed below. The file *opt.nio* is produced which contains 7 columns: energy (Ry), followed by the diagonal elements of the dielectric tensor (x,y,z components) for both spin channels. This can be compared with *EPS0001.dat*

Next we produce the files necessary for the **GW** calculation

```
$ echo -1 | lmfgwd `cat switches-for-lm` ctrl.nio
```

If we run **lmgwsc** as is then **lmgwsc** will perform a **GW** calculation with the **W** calculated within the RPA. To include ladder diagrams in **W** we need to edit *GWinput* and change/add the token UseBSEW 1

```
$ nano GWinput
UseBSEW 1
NumCondStates 11
NumValStates 16
```

The NumValStates and NumCondStates categories are the number of valence and conduction states considered in the BSE. Since the BSE is rather expensive (in compuational time and memory) all transitions between states cannot be considered at this level. States within a particular energy window about the Fermi level are considered in the BSE; the rest are treated at the level of the RPA. For this tutorial, 16 Valence and 11 conduction states are enough to demonstrate the method. There are other parameters that can be changed that affect the BSE program and some of these (along with their defaults and a description) are listed below:

```
ImagEner1 0.02 ! Broadening in Rydberg (~0.25eV)
ImagEner2 0.02 ! Broadening at Max Omega (see file freq_r). Broadening increases linearly from ImagEner1 at omega=0
NumValStates2 Nv2 ! Number of valence states in BSE for spin channel 2
NumCondStates2 Nc2 ! As above, but for the conduction states
Gauss_BSE off ! if on, use Gaussian broadening for the imaginary part of epsilon.
Gauss_cutoff 0.0 ! Gaussians used up to this cut-off, ensures the peaks aren't broadened too much near omega=0
```

Next we perform one iteration of QS**GW** with ladder diagrams now included in the screening.

```
$ lmgwsc --iter=1 --maxit=1 --code2 --openmp=N --sym --getsigp --insul=22 nio > out.lmgwsc_bse
```

The switches are the usual **lmgwsc** ones. **N** is the number of processors you wish to use in the calculation. In the file *out.lmgwsc_bse* we see that the output from the bethesalpeter code is directed to the file *lbse*. The RMS change in Sigma (see the bottom of *out.lmgwsc_bse*) is about 5.3E-3. This is a way of comparing H_{0} and H, i.e., comparing the QS**GW** Sigma without ladders in **W** and the case with ladders. Another iteration of **lmgwsc** reduces the calculated bandgap from 4.61 eV (see below) to 4.45 eV, and this is done by changing the switch **–maxit=2** in the command above.

Using the new sigm file created we can now run lmf and see how the band gap has changed.

```
$ lmf -vnit=12 --iactiv=no `cat switches-for-lm` ctrl.nio > out.lmf_QSGW_BSE
```

We find that the gap has reduced to 4.61 eV (from 4.86 eV). We can now produce a band structure from this converged density.

### Band Structure

The band structure can be produced in the usual way. In this example we follow the method on this page to produce a plot comparing band structures at different levels of theory (in this case QS**GW**(RPA) and BSE on top of QS**GW**). Note that we use ghostscript to view the postscript:

```
$ lmf -vnit=1 --rs=1,0 `cat switches-for-lm` --band:fn=syml ctrl.nio
$ cp bnds.nio bnds.bse
$ echo -9,16 / | plbnds -wscl=1,.8 -fplot -ef=0 -scl=13.6 --nocol -lbl=X,G,L,U,T -spin2 -dat=bse bnds.bse
$ echo -9,16 / | plbnds -wscl=1,.8 -fplot -ef=0 -scl=13.6 --nocol -dat=qsgw -spin2 bnds.qsgw
$ echo "% char0 colr=3,bold=5,clip=1,col=1,.2,.3" >>plot.2bands
$ echo "% char0 colb=2,bold=4,clip=1,col=.2,.3,1" >>plot.2bands
$ awk '{if ($1 == "-colsy") {sub("-qr","-lt {colg} -qr");sub("dat","green");sub("green","bse");sub("colg","colr");print;sub("qsgw","bse");sub("colr","colb");print} else {print}}' plot.plbnds >> plot.2bands
$ fplot -f plot.2bands
$ gs fplot.ps
```

The plot below has the vertex corrected QS**GW** in red (after one iteration of QSGW with BSE) and bare QS**GW** in blue.

### Optical absorption spectrum

See this page for a description of how to calculate an optical absorption spectrum for LiF. A list of commands are presented below for producing one for NiO. Edit the **OPTICS** catergory in *ctrl.nio* to calculate the transition dipole matrix elements (TDMEs) needed for the macroscopic dielectric function

```
$ nano ctrl.nio
OPTICS
MODE=1
LTET=0
WINDOW=0 4
NPTS=1001
FILBND=1 22
EMPBND=23 100
MEFAC=1
```

**MEFAC**=1 will include a contribution to the TDMEs from the self-energy (see Phys. Rev. Mat. 2, 034603). Note that **lmf** cannot run with **MEFAC**=1 and -vmetal=5 (in *switches-for-lm*) together. Therefore we change -vmetal to 3.

The set of commands to calculate the TDMEs are

```
$ echo 6 | lmfgwd `cat switches-for-lm` ctrl.nio # Ensures q-points and phases are the same between lmf and lmgw
$ lmf -vmetal=3 --opt:woptmc:permqp --revec:gw -vnit=1 --rs=1,0 -vloptic=1 `cat switches-for-lm` ctrl.nio
$ cp optdatac.nio optdatabse # optdatabse is the file read by bethesalpeter
```

Edit the *GWinput* for the calculation of the absorption:

```
$ nano GWinput
NumCondStates 11 !
Qvec_bse 0 0 1 ! q-vector direction
NumValStates 16
MaxOmega_BSE 2.0
MinOmega_BSE 0.0
numOmega_BSE 2001
Gauss_BSE on
ImagEner1 0.05 ! Large broadening to match experiment!
```

Now run

```
$ echo -1 | /path-to-questaal/lm/build/code2/bethesalpeter > out.rpa # Produces an RPA absorption spectrum, which can be compared with that in EPS0001.dat
$ cp eps_BSE.out eps_RPA.out
$ export MKL_NUM_THREADS=N # Diagonalization routine should run quicker
$ echo 0 | /path-to-questaal/lm/build/code2/bethesalpeter > out.bse
```

The absorption spectrum is written to the file *eps_BSE.out*, with 3 columns: energy (eV) followed by the real and imaginary part of the macroscopic dielectric function. Note the blue shift w.r.t. experiment, which can be reduced if we further iterate the QS**GW** including ladder diagrams, and if we could include the electron-phonon contribution to **W**. We plot eps_BSE.out against the experiment and the RPA QS**GW** spectrum from the file *EPS00001.dat*