# 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 $(G^0W^0)$ calculations are perturbations to a DFT calculation (such as LDA). They are simpler than QSGW calculations, because only the diagonal part of $\Sigma^0$ 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. $\Sigma^0$ 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 $\Gamma$ 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 $\Gamma$ 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 ($\Gamma$ 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 $\Gamma$ 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 $\Sigma$. 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 ($\Sigma$-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 $\Gamma$ (state 4) is -0.29 eV. The difference between the two (0.58 eV) is the LDA $\Gamma$-X bandgap. The quasi-particle energies including a Z factor are listed in the next column eQP. The quasipartcle $\Gamma$-X bandgap (1.08 eV) is given by the difference betwee state 5 in the X block of points and state 4 in the $\Gamma$ 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 $\Gamma$-X gap from the QSGW tutorial is around 1.4 eV. The QSGW approximation is known to systematically overestimates band gaps for most semiconductors.

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