# Introductory GW tutorial: LDA-based GW for Si

This tutorial begins with an LDA calculation for Si, starting from an init file. Following this is a demonstration of a 1-shot GW calculation. Click on the **GW Theory** dropdown menu below for a brief description of the 1-shot GW scheme. A complete summary of the commands used throughout is provided in the **Commands summary** dropdown menu. Theory for GW and QSGW, and its implementation in the Questaal suite, can be found in Phys. Rev. B76, 165106 (2007).

### Theory

1-shot GW calculations are perturbations to a DFT calculation (such as LDA). They are simpler than QSGW calculations, because only the diagonal part of is normally calculated (this is an approximation) and only one self-energy is calculated (single iteration). On the other hand, 1-shot calculations are sensitive to the starting point and you do not have the luxury of interpolating between k points to get full bandstructures, as is the case in QSGW. As a result, it is only possible to calculate 1-shot corrections for k points that lie on the k-point mesh used in the self-energy calculation.

The self-energy enters the Hamiltonian as a perturbation and gives us access to quasi-particle (QP) energies. The QP energies are the main output of a 1-shot GW calculation.

The DFT executable is **lmf**. **lmfgwd** is similar to **lmf**, but it is a driver whose purpose is to set up inputs for the *GW* code. is made by a shell script **lmgw1-shot**.

### Command summary

```
blm init.si --express --gmax=5 --nk=4 --nit=20 --gw #use blm tool to create actrl and site files
cp actrl.si ctrl.si && lmfa si && cp basp0.si basp.si #copy actrl to recognised ctrl prefix, run lmfa and copy basp
lmf ctrl.si > out.lmfsc #make self-consistent
echo -1 | lmfgwd si #make GWinput file
nano GWinput #edit k-point lines
lmgw1-shot --ht si #1-shot GW calculation
```

### LDA calculation

To carry out a self-consistent LDA calculation, we use the lmf code. The steps follow those from the lmf and QSGW tutorials; you should refer to these for additional details. Copy the following lines to a file called *init.si*:

```
LATTICE
ALAT=10.26
PLAT= 0.00000000 0.50000000 0.50000000
0.50000000 0.00000000 0.50000000
0.50000000 0.50000000 0.00000000
# pos means cartesian coordinates, units of alat
SITE
ATOM=Si POS= 0.00000000 0.00000000 0.00000000
ATOM=Si POS= 0.25000000 0.25000000 0.25000000
```

Run the following commands to obtain a self-consistent density:

```
blm init.si --express --gmax=5 --nk=4 --nit=20 --gw
cp actrl.si ctrl.si && lmfa si && cp basp0.si basp.si
lmf ctrl.si > out.lmfsc
```

Note that we have included an extra **--gw** switch, which tailors the *ctrl* file for a *GW* calculation. To see how it affects the *ctrl* file, try running blm without **--gw**. The switch modifies the basis set section (see the **autobas** line) to increase the size of the basis, which is necessary for *GW* calculations.

After running these commands, we now have a self-consistent LDA density. Check the *out.lmfsc* file and you should find a converged gap of around 0.58 eV. Now that we have the input eigenfunctions and eigenvalues, the next step is to carry out a *GW* calculation. For this, we need an input file for the *GW* code.

### Making GWinput

As in the QSGW case, a template *GWinput* file is made by running the following command:

```
echo -1 | lmfgwd si #make GWinput file
```

Take a look at *GWinput*, the k mesh is specified by n1n2n3 in the following line:

```
n1n2n3 4 4 4 ! for GW BZ mesh
```

In the QSGW tutorial we changed the mesh to 3x3x3 to speed things up. This time we will use the 4x4x4 mesh as the 3x3x3 mesh does not include the X and L points that are of particular interest. You may want to review the QSGW tutorial for a discussion of k-point convergence. The number of states (bands) to consider is specified in the following section:

```
*** no. states and list of band indices to make Sigma and QP energies
8
1 2 3 4 5 6 7 8
```

Just below the no. of states is a section that specifies the k points to consider:

```
*** q-points (must belong to mesh of points in BZ).
3
1 0.0000000000000000 0.0000000000000000 0.0000000000000000
2 -0.2500000000000000 0.2500000000000000 0.2500000000000000
3 -0.5000000000000000 0.5000000000000000 0.5000000000000000
4 0.0000000000000000 0.0000000000000000 0.5000000000000000
5 -0.2500000000000000 0.2500000000000000 0.7500000000000000
6 -0.5000000000000000 0.5000000000000000 1.0000000000000000
7 0.0000000000000000 0.0000000000000000 1.0000000000000000
8 0.0000000000000000 0.5000000000000000 1.0000000000000000
```

These are the 8 irreducible k points of the 4x4x4 mesh, including X (0,0,1) and L (-1/2,1/2,1/2). You can calculate QP corrections for all of these points but we will only calculate QP corrections at and X in this tutorial. The number **3** just below the q-points line tells the GW codes how many points to calculate QP corrections for. Use a text editor to change this number to **2** and move line 7 (which contains the X point) to appear second. The first few lines of your file should look like this:

```
*** q-points (must belong to mesh of points in BZ).
2
1 0.0000000000000000 0.0000000000000000 0.0000000000000000
7 0.0000000000000000 0.0000000000000000 1.0000000000000000
3 -0.5000000000000000 0.5000000000000000 0.5000000000000000
```

As specified above, the GW code will only calculate QP corrections for the first two lines (those containing and X).

### Running 1-shot GW

We can now run the 1-shot GW calculation, this is done with the **lmgw1-shot** script:

```
lmgw1-shot --ht si #1-shot GW calculation
```

The **insul** switch tells the code where to find the valence band maximum (8 electrons and spin degenerate so the fourth band is the valence band). Take a look at the **lmgw1-shot** output:

```
lmgw 16:22:08 : invoking /users/ms4/bin/lmfgwd --jobgw=0 --gwcode=2 --no-iactive si >llmfgw00
lmgw 16:22:08 : invoking /users/ms4/bin/code2/qg4gw --job=1 >lqg4gw
lmgw 16:22:08 : invoking /users/ms4/bin/lmfgwd --jobgw=1 --gwcode=2 --no-iactive si >llmfgw01
lmgw 16:22:09 : invoking /users/ms4/bin/lmf2gw_2 si >llmf2gw
lmgw 16:22:09 : invoking /users/ms4/bin/code2/rdata4gw_v2 --job=0 >lrdata4gw
lmgw 16:22:09 : invoking /users/ms4/bin/code2/heftet --job=1 >leftet
lmgw 16:22:09 : invoking /users/ms4/bin/code2/hchknw --job=0 >lchknw
lmgw 16:22:09 : invoking /users/ms4/bin/code2/hbasfp0 --job=3 >lbasC
lmgw 16:22:09 : invoking /users/ms4/bin/code2/hvccfp0 --job=0 >lvccC
lmgw 16:22:14 : invoking /users/ms4/bin/code2/hsfp0 --job=3 >lsxC
lmgw 16:22:15 : invoking /users/ms4/bin/code2/hbasfp0 --job=0 >lbas
lmgw 16:22:15 : invoking /users/ms4/bin/code2/hvccfp0 --job=0 >lvcc
lmgw 16:22:20 : invoking /users/ms4/bin/code2/hsfp0 --job=11 >lsx
lmgw 16:22:21 : invoking /users/ms4/bin/code2/hx0fp0 --job=11 >lx0
lmgw 16:22:30 : invoking /users/ms4/bin/code2/hsfp0 --job=12 >lsc
lmgw 16:22:38 : invoking /users/ms4/bin/code2/hqpe --job=0 >lqpe
```

The first few lines are preparatory steps, the main *GW* calculation begins on the line containing the file name **lbasC**. The three lines with **lbasC**, **lvccC** and **lsxC** are the steps that calculate the core contributions to the self-energy. The following lines up to the one with **lsc** are for the valence contribution to the self-energy. The **lsc** step, calculating the correlation part of the self-energy, is usually the most expensive part. The **lqpe** line assembles components of the potential and makes the **QPU** file. Further information can be found in the annotated GW output page.

The resulting quasi-particle (QP) energies are reported in the **QPU** file. Your **QPU** file should look like this (note that the last two columns have been removed in the output below):

```
q state SEx SExcore SEc vxc dSE dSEnoZ eLDA eQP eQPnoZ eHF Z
0.00000 0.00000 0.00000 1 -17.40 -1.81 6.70 -12.47 -0.02 -0.04 -12.27 -12.29 -12.31 -19.01 0.65
0.00000 0.00000 0.00000 2 -12.92 -1.96 1.30 -13.62 0.02 0.03 -0.29 -0.27 -0.26 -1.56 0.78
0.00000 0.00000 0.00000 3 -12.92 -1.96 1.30 -13.62 0.02 0.03 -0.29 -0.27 -0.26 -1.56 0.78
0.00000 0.00000 0.00000 4 -12.92 -1.96 1.30 -13.62 0.02 0.03 -0.29 -0.27 -0.26 -1.56 0.78
0.00000 0.00000 0.00000 5 -5.56 -1.42 -4.01 -11.82 0.63 0.82 2.22 2.85 3.04 7.05 0.77
0.00000 0.00000 0.00000 6 -5.56 -1.42 -4.01 -11.82 0.63 0.82 2.22 2.85 3.04 7.05 0.77
0.00000 0.00000 0.00000 7 -5.56 -1.42 -4.01 -11.82 0.63 0.82 2.22 2.85 3.04 7.05 0.77
0.00000 0.00000 0.00000 8 -5.85 -3.72 -4.57 -15.20 0.80 1.06 2.94 3.75 4.00 8.58 0.76
0.00000 0.00000 1.00000 1 -15.93 -2.11 4.80 -13.20 -0.03 -0.04 -8.13 -8.16 -8.17 -12.97 0.69
0.00000 0.00000 1.00000 2 -15.93 -2.11 4.80 -13.20 -0.03 -0.04 -8.13 -8.16 -8.17 -12.97 0.69
0.00000 0.00000 1.00000 3 -13.35 -1.69 2.30 -12.59 -0.12 -0.16 -3.16 -3.28 -3.32 -5.61 0.74
0.00000 0.00000 1.00000 4 -13.35 -1.69 2.30 -12.59 -0.12 -0.16 -3.16 -3.28 -3.32 -5.61 0.74
0.00000 0.00000 1.00000 5 -5.04 -0.92 -3.63 -10.25 0.52 0.66 0.29 0.81 0.95 4.58 0.79
0.00000 0.00000 1.00000 6 -5.04 -0.92 -3.63 -10.25 0.52 0.66 0.29 0.81 0.95 4.58 0.79
0.00000 0.00000 1.00000 7 -3.67 -2.35 -6.62 -13.50 0.63 0.86 9.73 10.35 10.59 17.20 0.73
0.00000 0.00000 1.00000 8 -3.67 -2.35 -6.62 -13.50 0.63 0.86 9.73 10.35 10.59 17.20 0.73
```

In the *GWinput* file we specified that the QP energies are to be calculated at two k points ( and X) and for 8 states each. The first three columns above list the k point x, y and z components. There is a block of 8 points, a space and then a block of 8 X points.

The **SEx**, **SExcore** and **SEc** columns contain the exchange and correlation terms, with the exchange divided into core and valence parts. These terms add to give you the GW self-energy . In the first line above, the three values sum to give a self-energy of around -12.51 eV. **vxc** is the matrix element of the LDA exchange-correlation potential for a given q-point and state, it is subtracted from the GW self-energy to get the QP shifts; you are adding (-vxc) to the LDA single-particle hamiltonian. Subtracting **vxc** from -12.51 gives you the -0.04 eV shift shown in the **dSEnoZ** column. This is the QP shift relative to LDA without a Z factor. The Z factor is a correction term that accounts for the fact that the self-energy is evaluated at the LDA eigenvalue (it should be evaluated at the QP eigenvalue but this is not known in advance). The **dSE** column is the QP shift relative to LDA including a z factor, it is obtained by multiplying the **dSEnoZ** value by the Z factor found in the **Z** column.

The column labelled **eLDA** contains the LDA eigenvalues. The conduction band energy (state 5) at the X point is 0.29 eV and the valence band energy at (state 4) is -0.29 eV. The difference between the two (0.58 eV) is the LDA -X bandgap. The quasi-particle energies including a Z factor are listed in the next column **eQP**. The quasipartcle -X bandgap (1.08 eV) is given by the difference betwee state 5 in the X block of points and state 4 in the block of points. The **eQPnoZ** column contains the QP energies without a Z factor correction. Lastly, the **eHF** column contains the Hartree-Fock eigenvalue energies which are obtained by subtracting the **vxc** value and adding the exchange terms **SEx** and **SExcore**.

The 1-shot bandgap is an improvement over the LDA, but it still underestimates the experimental value of 1.32 eV (at 0 K). This tendency is a common feature of semiconductors. It was not recognized for a long time because pseudopotential calculations (nearly all calculations were pseudopotential based until about 2005) tend to put levels too high, and somewhat remarkably, yielded fortutitously good agreement with experiment in many cases. On the other hand, the -X gap from the QSGW tutorial is around 1.4 eV. The QSGW approximation is known to systematically overestimates band gaps for most semiconductors.

### Additional Exercises

**1) GaAs**

Try a 1-shot GW calculation for GaAs. Note that the code automatically treats the Ga d state as valence (adds a local orbital). This requires a larger GMAX. You also need to run lmfa a second time to generate a starting density that includes this local orbital. The lmfa line for GaAs should be:

```
$ lmfa ctrl.gas; cp basp0.gas basp.gas; lmfa ctrl.gas
```

The init file is:

```
# init file for GaAs
LATTICE
# SPCGRP=
ALAT=10.69
PLAT= 0.00000000 0.50000000 0.50000000
0.50000000 0.00000000 0.50000000
0.50000000 0.50000000 0.00000000
# 2 atoms, 1 species
SITE
ATOM=Ga POS= 0.00000000 0.00000000 0.00000000
ATOM=As POS= 0.25000000 0.25000000 0.25000000
```