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).
Remark: This tutorial is written for the classic GW version of the code. It will be updated in due course.
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 --insul=4 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
See GWinput documentation for description of its structure and tags. Also, the QSGW tutorial has a bit more description of GWinput.
Take a look at GWinput; note in particular the k mesh is specified by tag n1n2n3 in the following line:
n1n2n3 4 4 4 ! for GW BZ mesh
The QSGW tutorial it suggests you change the mesh 3×3×3 to speed things up. It can do this because QSGW mode has the ability to interpolate the self-energy to any k point. The equivalent facility is not available in 1-shot mode.
This time we will use the 4×4×4 mesh as the 3×3×3 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 1but we will only calculate QP corrections at , L, 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 move line 7 (which contains the X point) to appear third. The first few lines of your file should look like this:
*** q-points (must belong to mesh of points in BZ).
3
1 0.0000000000000000 0.0000000000000000 0.0000000000000000
3 -0.5000000000000000 0.5000000000000000 0.5000000000000000
7 0.0000000000000000 0.0000000000000000 1.0000000000000000
As specified above, the GW code will only calculate QP corrections for the first three lines (those containing , L, 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 --insul=4 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=4 >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 often most expensive part, depending on how many states you ask for. 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 quasi-particle (QP) shifts and energies are reported in the QPU file; its structure is documented here. Your file should look like this:
q state SEx SExcore SEc vxc dSE dSEnoZ eLDA eQP eQPnoZ eHF Z FWHM=2Z*Simg ReS(elda)
0.00000 0.00000 0.00000 1 -17.40 -1.82 6.70 -12.47 -0.03 -0.04 -11.98 -12.02 -12.04 -19.01 0.65 1.25745 -12.51119
0.00000 0.00000 0.00000 2 -12.92 -1.96 1.30 -13.62 0.02 0.03 -0.00 -0.00 0.00 -1.56 0.78 0.00000 -13.59116
0.00000 0.00000 0.00000 3 -12.92 -1.96 1.30 -13.62 0.02 0.03 -0.00 0.00 0.00 -1.56 0.78 0.00000 -13.59117
0.00000 0.00000 0.00000 4 -12.92 -1.96 1.30 -13.62 0.02 0.03 0.00 0.00 0.00 -1.56 0.78 0.00000 -13.59117
0.00000 0.00000 0.00000 5 -5.56 -1.42 -4.01 -11.82 0.63 0.82 2.51 3.12 3.30 7.05 0.77 -0.00096 -10.99915
0.00000 0.00000 0.00000 6 -5.56 -1.42 -4.01 -11.82 0.63 0.82 2.51 3.12 3.30 7.05 0.77 -0.00096 -10.99916
0.00000 0.00000 0.00000 7 -5.56 -1.42 -4.01 -11.82 0.63 0.82 2.51 3.12 3.30 7.05 0.77 -0.00096 -10.99916
0.00000 0.00000 0.00000 8 -5.85 -3.72 -4.57 -15.20 0.80 1.06 3.23 4.01 4.26 8.57 0.76 -0.00274 -14.13973
-0.50000 0.50000 0.50000 1 -16.79 -2.09 5.68 -13.15 -0.04 -0.06 -9.64 -9.70 -9.72 -15.67 0.72 0.98211 -13.20291
-0.50000 0.50000 0.50000 2 -14.81 -1.66 4.27 -12.12 -0.06 -0.08 -7.02 -7.10 -7.12 -11.66 0.71 0.39708 -12.19903
-0.50000 0.50000 0.50000 3 -13.14 -1.92 1.70 -13.30 -0.05 -0.06 -1.21 -1.27 -1.29 -3.25 0.76 0.00000 -13.36701
-0.50000 0.50000 0.50000 4 -13.14 -1.92 1.70 -13.30 -0.05 -0.06 -1.21 -1.27 -1.29 -3.25 0.76 0.00000 -13.36702
-0.50000 0.50000 0.50000 5 -5.81 -2.16 -3.86 -12.66 0.65 0.83 1.42 2.05 2.23 5.82 0.78 -0.00000 -11.83023
-0.50000 0.50000 0.50000 6 -4.90 -0.99 -4.26 -10.98 0.63 0.83 3.28 3.90 4.08 8.08 0.77 -0.00729 -10.15165
-0.50000 0.50000 0.50000 7 -4.90 -0.99 -4.26 -10.98 0.63 0.83 3.28 3.90 4.08 8.08 0.77 -0.00729 -10.15166
-0.50000 0.50000 0.50000 8 -2.38 -0.59 -5.33 -8.87 0.44 0.57 7.58 8.00 8.12 13.19 0.77 -0.32644 -8.29502
0.00000 0.00000 1.00000 1 -15.93 -2.12 4.80 -13.20 -0.03 -0.04 -7.84 -7.89 -7.91 -12.97 0.69 0.54541 -13.24475
0.00000 0.00000 1.00000 2 -15.93 -2.12 4.80 -13.20 -0.03 -0.04 -7.84 -7.89 -7.91 -12.97 0.69 0.54540 -13.24471
0.00000 0.00000 1.00000 3 -13.35 -1.69 2.30 -12.59 -0.12 -0.16 -2.87 -3.01 -3.05 -5.61 0.74 0.00361 -12.74128
0.00000 0.00000 1.00000 4 -13.35 -1.69 2.30 -12.59 -0.12 -0.16 -2.87 -3.01 -3.05 -5.61 0.74 0.00361 -12.74128
0.00000 0.00000 1.00000 5 -5.04 -0.92 -3.63 -10.25 0.52 0.66 0.58 1.08 1.21 4.58 0.79 -0.00000 -9.59027
0.00000 0.00000 1.00000 6 -5.04 -0.92 -3.63 -10.25 0.52 0.66 0.58 1.08 1.21 4.58 0.79 -0.00000 -9.59034
0.00000 0.00000 1.00000 7 -3.67 -2.35 -6.62 -13.50 0.63 0.86 10.02 10.62 10.85 17.20 0.73 -0.64702 -12.64267
0.00000 0.00000 1.00000 8 -3.67 -2.35 -6.62 -13.50 0.63 0.86 10.02 10.62 10.85 17.20 0.73 -0.64702 -12.64267
In the GWinput file we specified that the QP energies are to be calculated at three k points (, L, 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. However, there is a good reason to omitt the Z factor; it is a kind of poor man’s self-consistency. See, for example the Appendix in Phys. Rev. B74, 245125.
The column labelled eLDA contains the LDA eigenvalues. The conduction band energy (state 5) at the X point is 0.58 eV and the valence band energy at (state 4) is 0 eV. (It is zero because command-line argument --insul=4 told lmgw to use the fourth band as the energy zero. 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