Introductory QSGW 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 (QSGW) calculation. An example of the 1-shot GW code is provided in a separate tutorial. Click on the dropdown menu below for a brief description of the QSGW scheme. A complete summary of the commands used throughout is provided in a separate dropdown menu. Theory for GW and QSGW, and its implementation in the Questaal suite, can be found in Phys. Rev. B76, 165106 (2007).


QSGW flowchart

Each iteration of a QSGW 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 quaisparticle self-consistency, the GW code makes a “quasiparticlized” self-energy Σ0, which is used to construct the effective one-body hamiltonian for the next cycle. The process is iterated until the change in Σ0 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. Σ0 is made by a shell script lmgw. The entire cycle is managed by a shell script lmgwsc.

Before any self-energy Σ0 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 G[Σ(ω)] and non-interacting G00]. At self-consistency the poles of G and G0 coincide: this is a unique and very advantageous feature of QSGW. 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 QSGW 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 Σ0. This is called the “0th iteration.” If only the diagonal parts of Σ0 are kept, the “0th” iteration corresponds to what is sometimes called 1-shot GW, or as GLDAWLDA.

In the next iteration, Σ0VxcLDA 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 VxcLDA cancels the exchange-correlation potential from the LDA hamiltonian.) This process is repeated until the RMS change in Σ0 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.


Command summary

$ 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
$ vim GWinput                                              #change GW k mesh to 3x3x3
$ 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=3 && 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. One thing to note is the NKABC= token, which defines the GW k-point mesh. It is specified in the same way as the lower case nkabc for the LDA calculation.

Now check the output file out.lmfsc. The self-consistent gap is reported to be around 0.58 eV as can be seen by searching for the last occurence of the word ‘gap’. Note that this result differs slightly to 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 QSGW) uses one main user-supplied input file, GWinput. The script lmfgwd can create a template GWinput file for you by running the following command:

$ echo -1 | lmfgwd si                              #make GWinput file

The lmfgwd script 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’. Take a look at GWinput, it is a rather complicated input file but we will only consider the GW k-point mesh for now (further information can be found on the GWinput page). The k mesh is specified by n1n2n3 in the GWinput file, look for the following line:

$ n1n2n3  4  4  4 ! for GW BZ mesh

When creating the GWinput file, lmfgwd checks the GW section of the ctrl file for default values. The ‘NKABC= 4’ part of the DFT input file (ctrl.si) is read by lmfgwd and used for n1n2n3 in the GW input file. Remember if only one number is supplied in NKABC then that number is used as the division in each direction of the reciprocal lattice vectors, so 4 alone means a 4 x 4 x 4 k mesh. To make things run 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. As is the case with the LDA, it is very important to control k convergence. However, a coarser mesh can often be used in GW because the self-energy generally varies much more smoothly with k than does the kinetic energy. This is fortunate because GW calculations are much more expensive. It is important to note that convergence tests will have to be performed for any new system. These can be time consuming and unfortunately there are no shortcuts.


Running QSGW


We are now ready for a QSGW calculation, this is run using the shell script lmgwsc:

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

The 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 QSGW 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 QSGW 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 self-energy is converted into an effective exchange-correlation potential in the lqpe step and the final few lines are to do with its handling. Further information can be found in the annotated GW output page.

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 calculation but this the zeroth iteration GW potential is used. The following line in the llmf file 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 RDSIG variable in the ctrl file. 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.

lmgwsc : iter 3 of 5  RMS change in sigma = 1.19E-05  Tolerance = 2e-5

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

Inspect the lmf output and you can find that the gap is now around 1.42 eV. The --rs option tells lmf to read from the restart file, which contains the 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.

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.


Additional Exercises

1) Correct gap

This is actually the Γ-X gap; the true gap is around 1.30 eV 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 4 x 4 x 4 k mesh.

3) Plotting bands

Plotting bands is the same as in the DFT case. Try repeating the same steps for your QSGW converged density.

4) GaAs

Try a QSGW 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

Edit This Page