# Introductory QS*GW* Tutorial

This tutorial begins with an LDA calculation for Si, starting from an *init* file. Following this is a demonstration of a quasi-particle self-consistent GW (QS*GW*) calculation. An example of the 1-shot GW code is provided in a separate tutorial.

After the potential is obtained, the energy bands are drawn, and the macroscopic dielectric function is calculated in the RPA. This tutorial explains how to calculate with ladder diagrams.

Click on the Command Summary below for a summary of the commands explained in this tutorial. Theory for *GW* and QS*GW*, and its implementation in the Questaal suite, can be found in Phys. Rev. B76, 165106 (2007).

Each iteration of a QS*GW* calculation has two main parts: a section that uses effective one-body hamiltonians to make the density *n* (as in DFT), and the *GW* code that makes the self-energy of an interacting hamiltonian. For quasiparticle self-consistency, the *GW* code makes a “quasiparticlized” self-energy , which is used to construct the effective one-body hamiltonian for the next cycle. The process is iterated until the change in becomes small.

The one-body executable is **lmf**. The script **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 **lmgw**. The entire cycle is managed by a shell script **lmgwsc**, which runs the cycle shown in the figure above, and also does some bookkeeping, such as monitoring connvergence to self-consistency.

Before any self-energy exists, it is assumed to be zero. Thus the one-body hamiltonian is usually the LDA, though it can be something else, e.g. LDA+U.

*Note:* in some circumstances, e.g. to break time reversal symmetry inherent in the LDA, you need to start with LDA+U.

Thus, there are two self-energies and two corresponding Green’s functions: the interacting and non-interacting . At self-consistency the poles of and coincide: this is a unique and very advantageous feature of QS*GW*. It means that there is no “mass renormalization” of the bandwidth — at least at the *GW* level.

Usually the interacting isn’t made explicitly, but you can do so, as explained in this tutorial.

In short, a QS*GW* calculation consists of the following steps. The starting point is a self-consistent DFT calculation (usually LDA). The DFT eigenfunctions and eigenvalues are used by the GW code to construct a self-energy . This is called the “0^{th} iteration.” If only the diagonal parts of are kept, the “0^{th}” iteration corresponds to what is sometimes called 1-shot *GW*, or as *G*^{LDA}W^{LDA}.

In the next iteration, is added to the LDA hamiltonian. The density is made self-consistent, and control is handed over to the *GW* part. (Note that for a fixed density cancels the exchange-correlation potential from the LDA hamiltonian.) This process is repeated until the RMS change in falls below a certain tolerance value. The final self-energy (QSGW potential) can be thought of as an effective exchange-correlation functional that is tailored to the system. This is very convenient as it can now be used in an analogous way to standard DFT to calculate properties such as the band structure.

```
nano init.si #create init file using lines from box below
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 si > out.lmfsc #make self-consistent
echo -1 | lmfgwd si #make GWinput file
lmgwsc --wt --insul=4 --tol=2e-5 --maxit=5 si #self-consistent GW calculation
vim ctrl.si #change number of iterations to 1
lmf si --rs=1,0 #lmf with QSGW potential to get QSGW band gap
lmgwclear #clean up directory
```

Alternatively, you can run the following two one-liners to get the same result. This assumes you have already created the *init* file.

```
blm init.si --express=0 --gmax=5 --nk=4 --nit=20 --gw --nkgw=4 && cp actrl.si ctrl.si && lmfa si && cp basp0.si basp.si && lmf si > out.lmfsc
echo -1 | lmfgwd si && lmgwsc --wt --insul=4 --tol=2e-5 --maxit=5 si > out.gwsc && lmgwclear && lmf si --rs=1,0 -vnit=1 > out.lmf_gwsc
```

### LDA calculation

The starting point is a self-consistent LDA density, you may want to review the DFT tutorial for silicon. 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 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.

Two new blocks of text, the **HAM** and **GW** categories, are also added towards the end of the *ctrl* file. The extra parameters in the **HAM** category handle the inclusion of a self-energy, actually a GW potential (see theory notes above), in a DFT calculation. The **GW** category provides some default values for parameters that are required in the *GW* calculation. The *GW* code has its own input file and the DFT *ctrl* file influences what defaults are set in it, we will come back to this later.

Now check the output file *out.lmfsc*. It should converge in 10 iterations as can be seen from the end of the file. Note that **ehf** and **ehk** should be essentially the same (-1156.140684 Ry).

Search for the last occurence of the word ‘**gap**’. It should be around 0.58. Note that this result differs slightly from that from the LDA tutorial because the **gw** switch increases the size of the basis set.

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

The *GW* package (both one-shot and QS*GW*) uses one main user-supplied input file, *GWinput*. The structure of *GWinput* is documented here.

**lmfgwd** supplies the interface connecting the LDA code to the *GW* package. (The *GW* package only makes matrix elements of response functions or self-energies suitable for later analysis.)

It can also create a template *GWinput* file by running the following command:

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

This is particularly convenient because *GWinput* is not easy to read, or to construct by hand.

*Caution:* the *GWinput* file made by **lmfgwd** should be regarded only as template. It is up to you to ensure that the contents are reasonable, even while for the most the template is reasonable.

**lmfgwd** has multiple options and is designed to run interactively. Using ‘**echo -1**’ automatically passes it the ‘**-1**’ option that specifies making a template input file. You can try running it interactively by just using the command ‘**lmfgwd si**’ and then entering ‘**-1**’. Some tags in **GWinput** read from the **GW** category. To see which ones, try

```
lmfgwd --input | grep GW_
```

Tags in the **GW** category of *ctrl.si* are used only to modify *GWinput*; the *GW* code does not read from the *ctrl* file. All of the tags in the You don’t need to supply any tags; **lmfgwd** will use reasonable defaults. Take a look at *GWinput*; it is a rather complicated input file but most of the time it the template **lmfgwd** uses reasonable defaults.

The

Two tags you particularly need to monitor are **QpGcut_cou** and **n1n2n3**. The former controls the plane wave cutoff. The computational cost is senstive to this value, so you want to make it as low as possible without messing up the results. It is a good idea to try a couple of values to satisfy your self that you have made a sensible choice: make it as small as you can while keeping the QP levels reasonable well converged.

**n1n2n3** specifies the *GW* k-point mesh for now (further information can be found on the GWinput page). In *GWinput* look for the following line:

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

This number was culled from *ctrl.si*. When reading this file, **lmfgwd** looks for tag **GW_NKABC=4**. In this case only one one number is supplied; that number is used for all three of the reciprocal lattice vectors, so 4 alone implies a 4 x 4 x 4 k mesh. If you want to run the calculation a bit quicker, change the k mesh to 3 x 3 x 3 by editing the *GWinput* file line:

```
$ n1n2n3 3 3 3 ! for GW BZ mesh
```

The k mesh of 3 x 3 x 3 divisions is rough, but it makes the calculation fast and for Si the results are reasonable. In this tutorial we will stick with 4 divisions. In any case, you need to monitor the *k* convergence. These can be time consuming and unfortunately there are no shortcuts. Note however, that the *GW* mesh is distinct from the *k* mesh **lmf** uses (read from **BZ_NKABC**). Typically a coarser mesh can often be used for the *GW* self-energy because generally varies much more smoothly with *k* than does the kinetic energy. This is fortunate because GW calculations are much more expensive, and in the existing implementation the computation time scales as *n*^{2}.

Other convergence parameters to particularly watch out for are the product basis tolerance, the augmention lcutmx, and whether product basis functions are included for shallow core states.

### Running QSGW

We are now ready for a QS*GW* calculation, which you do with the shell script **lmgwsc**:

```
$ lmgwsc --wt --insul=4 --tol=2e-5 --maxit=0 si #zeroth iteration of QSGW calculation
```

Switch ‘**–wt**’ includes additional timing information in the printed output, **insul** refers to the number of occupied bands (normally spin degenerate so half the number of electrons), **tol** is the tolerance for the RMS change in the self-energy between iterations and **maxit** is the maximum number of QS*GW* iterations. Note that **maxit** is zero, this specifies that a single iteration is to be carried out starting from DFT with no self-energy (zeroth iteration).

Take a look at the line containing the file name *llmf*:

```
lmgw 15:26:47 : invoking mpix -np=8 /h/ms4/bin/lmf-MPIK --no-iactive cspi >llmf
```

Each QS*GW* iteration begins with a self-consistent DFT calculation by calling the program **lmf** and writing the output to the file *llmf*. We are starting from a self-consitent LDA density (we already ran **lmf** above) so this step is not actually necessary here. The next few lines are preparatory steps. The main *GW* calculation begins on the line containing the file name ‘**lbasC**’:

```
lmgw 16:27:55 : invoking /h/ms4/bin/code2/hbasfp0 --job=3 >lbasC
lmgw 16:27:55 : invoking /h/ms4/bin/code2/hvccfp0 --job=0 >lvccC ... 0.0m (0.0h)
lmgw 16:27:58 : invoking /h/ms4/bin/code2/hsfp0_sc --job=3 >lsxC ... 0.0m (0.0h)
lmgw 16:27:59 : invoking /h/ms4/bin/code2/hbasfp0 --job=0 >lbas
lmgw 16:27:59 : invoking /h/ms4/bin/code2/hvccfp0 --job=0 >lvcc ... 0.0m (0.0h)
lmgw 16:28:02 : invoking /h/ms4/bin/code2/hsfp0_sc --job=1 >lsx ... 0.0m (0.0h)
lmgw 16:28:02 : invoking /h/ms4/bin/code2/hx0fp0_sc --job=11 >lx0 ... 0.1m (0.0h)
lmgw 16:28:07 : invoking /h/ms4/bin/code2/hsfp0_sc --job=2 >lsc ... 0.1m (0.0h)
lmgw 16:28:13 : invoking /h/ms4/bin/code2/hqpe_sc 4 >lqpe
```

The three lines with **lbasC**, **lvccC** and **lsxC** are the steps that calculate the core contributions to the self-energy and 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 step.

The last step (**hpqe_sc**) collects terms to make quasiparticlized self-energy and writes it to file *sigm* for every irreducible **k** point. Actually it writes . This makes running **lmf** very convenient, since **lmf** simply has to add this term to the LDA potential. Further information can be found in the annotated GW output page.

The self-energy produced so far is essentially the same as *G*^{LDA}*W*^{LDA} generated in the one-shot tutorial, only now there is no Z factor and the full is generated. This is distinct from computing the level shift in 1first order perturbation theory as **lmgw1-shot** did (and most *GW* codes do). This requires only the diagonal , which is fast and easier to make (and is why QS*GW* is more expensive to do).

It is interesting to compare one-shot results from the 0th iteration to the output of **lmgw1-shot**.

```
$ lmf si
```

The following line in the standard output specifies that the *GW* potential is being used:

```
RDSIGM: read file sigm and create COMPLEX sigma(R) by FT ...
```

The *GW* potential is contained in the file *sigm*, **lmgwsc** also makes a soft link *sigm.si* so **lmf** can read it. The *GW* potential is automatically used if present, this is specified by the **HAM_RDSIG** tag in the *ctrl* file.

After the first band pass, **lmf** yields a gap of 1.21 eV, essentially identical to what **lmgw1-shot** gives without a Z factor. In general this is not true; but Si is very simple and the and LDA eigenfunctions are very similar.

Nevertheless, note that as the density is updated (the off-diagonal elements of mean that the eigenfunctions change), the gap increases to 1.26 eV. This can be expressed as a change , where is implicitly given from DFT through self-consistency in **lmf**.

Run the command again but this time set the number of iterations (**maxit**) to something like 5:

```
$ lmgwsc --wt --insul=4 --tol=2e-5 --maxit=5 si #self-consistent GW calculation
```

The iteration count starts from 1 since we are now starting with a self-energy from the zeroth iteration. Again, the iteration starts with a self-consistent DFT-like calculation, but now from zeroth iteration is added. Take a look at the *GW* output again and you can see that the rest of the steps are the same as before. After 3 iterations the RMS change in the self-energy is below the tolerance - the calculation is converged.

```
mgwsc : iter 3 of 5 RMS change in sigma = 5.14E-06 Tolerance = 2e-5 more=F Mon 21 May 2018 19:08:21 BST elapsed wall time 5.0m (0.1h) mark
```

Now that we have a converged self-energy (*sigm*) we can go back to using **lmf** to calculate additional properties. We only want to run a single iteration so change the number of iterations (nit) to 1 in the *ctrl* file. Run the following command:

```
$ lmf si --rs=1,0 #lmf with QSGW potential to get QSGW band gap
```

The `--rs`

switch tells **lmf** to read from the restart file, which contains the LDA density, but not to write to it (once we have a converged density we want to keep this fixed). More information on command line switches can be found here.

Inspect the **lmf** output and you can find that the gap is now around 1.33 eV. It is larger because the one-body hamiltonian generating has a wider gap, which which increases $W$ and thus .

Check your directory and you will see that a large number of files were created. The following command removes many redundant files:

```
$ lmgwclear #clean up directory
```

Further details can be found in the Additional exercises below.

### QSGW energy bands

**lmf** has a very powerful feature, that it can takes the inverse Bloch transform to put in real space, from the mesh of points it is calculated on From real space it can interpolate to any by performing a forward Bloch transform. The details are rather complicated, but they are explained in some detail in Section IIG of PRB76, 165106.

This feature enables us to compute the energy bands and any **k**, and allows us to draw the energy bands with minimal effort.

```
$ lmchk ctrl.si --syml~n=41~q=-.5,.5,.5,0,0,0,0,0,1
$ lmf ctrl.si --band~fn=syml
$ echo -13,10,5,10 | plbnds -fplot -ef=0 -scl=13.6 -lbl=L,G,X bnds.si
$ fplot -f plot.plbnds
```

**lmchk** makes a symmetry lines file, along the lines L-Γ and Γ-X. **lmf** generates the energy bands in symmetry-lines mode. Finally the last two commands convert *bnds.si* generated by **lmf** into a postscript figure.

### QSGW RPA dielectric function

Here we show how to calculate the macroscopic dielectric function for Si in the RPA. Another tutorial shows how to calculate the macroscopic dielectric function with ladder diagrams.

The RPA dielectric function is computed from a starting hamiltonian . You can use any available , e.g. LDA, Here we use the QS*GW* generated in the previous section.

To make a good dielectric function it is necesssary to modify *GWinput* a little.

The 4×4×4

*k*mesh used for QS*GW*is reasonably sufficient for*k*integrated propertiees, but it is not not sufficient to pick up the detailed structure of the macroscopic dielectric function. Edit the**n1n2n3**tag and increase the*k*mesh to 16×16×16 divisions:`n1n2n3 16 16 16 ! for GW BZ mesh`

For a good frequency resolution, shrink the frequency mesh spacing to

**dw**to 0.001 Hartree (or smaller, if desired)`dw 0.001 ! mesh spacing along Real axis (Hartree)`

The macroscopic dielectric function is . Because it takes an average over the unit cell, the detailed shape of inside the unit cell is not so important. You can reduce the cutoffs significantly with minimal change in . Reducing it to 2 makes negligible difference in this case:

`lcutmx(atom) = l-cutoff for the product basis 2 2`

Supply

*GWinput*with the*q*points for which you want to calculate . The code does not have the ability to take the limit analytically; you must choose a finite but small**q**. At the end of*GWinput*there are the following lines`QforEPSIBZ off <QforEPS> 0d0 0d0 0.015d0 </QforEPS>`

It tells the *GW* code to make just a one point, close to the point. This corresponds approximately to . If you want to calculate also, add a line **0.015d0 0d0 0d0**. Si is cubic and , so it is not necessary here.

You can either calculate the proper macroscopic dielectric function, , or one that neglects local fields, .

The latter is must faster, but it is approximate. Often the difference between the two is not large, however. Do one of:

```
lmgw --eps --ht si
lmgw --epsNLF --ht si
```

The latter writes a file *EPS0001.nlfc.dat*, while the former writes both that file and *EPS0001.dat*.

Th first few lines should like the following

```
q(1:3) w(Ry) eps epsi --- LFC included.
0.00000000 0.00000000 0.01500000 0.0000D+00 0.928993543761226D+01 0.726517122804292D-15 0.107643374565477D+00 -0.841822371139706D-17
0.00000000 0.00000000 0.01500000 0.1005D-02 0.928999688792784D+01 0.721333610616975D-15 0.107642662539476D+00 -0.835805128491751D-17
0.00000000 0.00000000 0.01500000 0.3025D-02 0.929049220411353D+01 0.727839656503130D-15 0.107636923645147D+00 -0.843253724471668D-17
0.00000000 0.00000000 0.01500000 0.5065D-02 0.929149658090234D+01 0.730784249155761D-15 0.107625288487475D+00 -0.846482210402463D-17
```

Colums 1-4 are the q-point and frequency, resepectively. Columns 5-6 are the real and imaginary parts of ; columns 7-8 are the real and imaginary parts of .

Note that even while Si has a fundamental gap of approximately 1.2 eV, it is an indirect gap semiconductor. Since the electron-phonon interaction is missing absorption from indirect transitions is not included and is zero until the first direct gap transition at 0.2516D Ry, or 3.42 eV.

You can read the static dielectric constant, from *EPS0001.dat*. The calculated value is 9.3, which is 79% of the actual (11.7 at 0K). This error largely originates from the neglect of ladder diagrams.

#### The AC conductivity

The conductivity and dielectric function are related as

What has been calculated above is the dielectric function relative to the vacuum permittivity, . Taking the real part of and scaling to appropriate units,

Using

we can write

Then

### Additional Exercises

**1) Fundamental gap**

The gap found in the tutorial is is actually the Γ-X gap; the true gap is a little smaller as can be seen by running **lmf** with a fine k mesh.

**2) Changing k-point mesh**

Test the convergence with respect to the *GW* *k* mesh by increasing to a 8 × 8 × 8 k mesh.

**3) GaAs**

Try a QS*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
```