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.

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:

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:

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:

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):

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: