# 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.

After the potential is obtained, the energy bands are drawn, and the macroscopic dielectric function $\varepsilon_M(\omega) = \left[\varepsilon^{-1}_{G=0,G^\prime=0}(q{\rightarrow}0,\omega)\right]^{-1}$ is calculated in the RPA. This tutorial explains how to calculate $\varepsilon_M(\omega)$ with ladder diagrams.

Theory for GW and QSGW, and its implementation in the Questaal suite, can be found in Phys. Rev. B76, 165106 (2007).

Self-consistent LDA calculation, followed by QSGW

nano init.si
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

echo -1 | lmfgwd si
lmgwsc --wt --insul=4 --tol=2e-5 --maxit=5 si
lmf si --rs=1,0
lmgwclear


QSGW energy bands (see tutorial to plot QSGW and LDA energy bands together)

lmchk ctrl.si  --syml~n=41~lblq:G=0,0,0,L=1/2,1/2,1/2,X=1,0,0,W=1,1/2,0,K=3/4,3/4,0~lbl=LGXWGK
lmf ctrl.si --quit=rho
lmf ctrl.si --band~fn=syml
echo -13,10,5,10 | plbnds -fplot -ef=0 -scl=13.6 bnds.si
fplot -f plot.plbnds
open fplot.ps


QSGW RPA dielectric function

lmgw --eps --ht si
mcx -f1p8e20.12 -r:s=1 EPS0001.dat -coll 4,5,6 -p -coll 1 -s13.6 -tog -coll 2:nc -ccat -inc 'x1>0&x1<10' > dat.lfc
fplot dat.lfc  -coll 3 -lt 3,bold=4 dat.lfc
open fplot.ps 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 $\Sigma(\omega)$ of an interacting hamiltonian. For quasiparticle self-consistency, the GW code makes a “quasiparticlized” self-energy $\Sigma^0$, which is used to construct the effective one-body hamiltonian for the next cycle. The process is iterated until the change in $\Sigma^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. $\Sigma^0$ 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 convergence to self-consistency.

This tutorial also uses some auxiliary executables, e.g. plbnds and fplot utilities. Postscript files are assumed to be rendered by the Mac-style open utility.

Before any self-energy $\Sigma^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. 1

Thus, there are two self-energies and two corresponding Green’s functions: the interacting $G[\Sigma(\omega)]$ and non-interacting $G^0[\Sigma^0]$. At self-consistency the poles of $G$ and $G^0$ 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 $\Sigma(\omega)$ 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 $\Sigma^0$. This is called the “0th iteration.” If only the diagonal parts of $\Sigma^0$ are kept, the “0th” iteration corresponds to what is sometimes called 1-shot GW, or as GLDAWLDA.

In the next iteration, $\Sigma^0-V_{xc}^\text{LDA}$ 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 $V_{xc}^\text{LDA}$ cancels the exchange-correlation potential from the LDA hamiltonian.) This process is repeated until the RMS change in $\Sigma^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.

### 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 occurrence 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. Before moving on to the GW calculation, let’s make the LDA energy bands. lmf requires a symmetry lines file. We’ll chose the popular sequence L,Γ,X,W,Γ,K. Make the file this way: $ lmchk ctrl.si  --syml~n=41~lblq:G=0,0,0,L=1/2,1/2,1/2,X=1,0,0,W=1,1/2,0,K=3/4,3/4,0~lbl=LGXWGK


Aside: If you need to see or remember what a particular command-line switch does, there is a web page that contains them all with descriptions organized by executable. For the –syml switch the command-line documentation just points to another page where it is explained in more detail.

Make the bands as follows:

$lmf ctrl.si --quit=rho # makes a pass without changing the density: sets Ef in moms.si$ lmf ctrl.si --band~fn=syml
$cp bnds.si bnds.blue # save it so we can compare to GW  Instead of plotting the LDA bands, we’ll defer plots until we have GW results to compare with. ### Making GWinput The GW package (both one-shot and QSGW) 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: $ lmfgwd --job=-1 si                              #make GWinput file


This is particularly convenient because GWinput is not easy to read, or to construct by hand. Its contents are documented on this page.

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 sensitive 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 n2.

Other convergence parameters to particularly watch out for are the product basis tolerance, the augmentation lcutmx, and whether product basis functions are included for shallow core states. The GWinput file is documented on this page.

Before running the QSGW, it’s a good idea to use lmfgwd to perform a sanity check on GWinput. It checks for the reasonableness of a number of of parameters, and whether sites in the same class have the same basis. Run it with job = -2:

$lmfgwd --job=-2 si #sanity check  ### Running QSGW We are now ready for a QSGW 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 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-consistent 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 $\Sigma^0$ and writes it to file sigm for every irreducible k point. Actually it writes $\Sigma^0-V_\mathrm{xc}^\mathrm{LDA}$. 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 GLDAWLDA generated in the one-shot tutorial, only now there is no Z factor and the full $\Sigma^{nn^\prime}$ 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 $\Sigma^{nn}$, which is fast and easier to make (and is why QSGW 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 --quit=rho  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. Before continuing, let’s make and preserve the energy bands from the 0th iteration. $ lmf ctrl.si --quit=rho                       # makes a pass without changing the density: sets Ef in moms.si
$lmf ctrl.si --band~fn=syml$ cp bnds.si bnds.red                          # save it so we can compare to other bands


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 $GW$ and LDA eigenfunctions are very similar. Another way to see this is to look how much the output density has changed. This line

 mixrho:  sought 2 iter from file mixm; read 0.  RMS DQ=7.44e-4


shows how much the density changes from the hamiltonian $H^\mathrm{LDA}{\to}H^\mathrm{LDA}+\Sigma-V^\mathrm{LDA}_\mathrm{xc}$. It is not a large change; this is to be expected because the LDA and GW eigenfunctions are similar, as they typically are with simple sp semiconductors.

The change is small but not negligible. Run lmf to self-consistency keeping the $\Sigma$ fixed, but allowing $V^\mathrm{LDA}$ to change with the density.

$rm mixm.si # remove the mixing file so as not to confuse lmf$ lmf si                                       # Update the density


Note that as the density is updated (the off-diagonal elements of $\Sigma^{nn^\prime}$ mean that the eigenfunctions change), the gap increases to 1.26 eV. This can be expressed as a change $(\delta V/\delta n) \Delta n = \chi^{-1} \Delta n$, where $\chi^{-1}$ is implicitly given from DFT through self-consistency in the lmf cycle.

#### Iterating QSGW to self-consistency

Run lmgwsc 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 $\Sigma^0-V_{xc}^\text{LDA}$ 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. lmgwsc : 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)  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 --quit=rho                                    #lmf with QSGW potential to get QSGW band gap


Since the change in Σ between iteration 2 to 3 is small, we expect the change in density to be small also. Indeed this is the case. From the lmf you should see something like

 mixrho:  sought 2 iter from file mixm; read 2.  RMS DQ=5.72e-6


Thus at self-consistency, Σ and the density both stop changing, ameliorating ambiguity in, e.g. the Fermi level. QSGW is not a conserving approximation, but it is fairly close to a conserving one, and indeed vastly better than one-shot GW. This points out one very useful ground state property QSGW can give you: the density. It can also give you spin densities, e.g. magnetic moments.

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 $\Sigma$ has a wider gap, which which increases $W$ and thus $\Sigma$.

Make the bands a third time so we can compare LDA, GW(iteration 0), and QSGW.

$lmf ctrl.si --quit=rho # makes a pass without changing the density: sets Ef in moms.si$ lmf ctrl.si --band~fn=syml
$cp bnds.si bnds.green # save it so we can compare to other bands  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 $\Sigma^0(\mathbf{q}$ to put $\Sigma^0$ in real space, from the mesh of points it is calculated on From real space it can interpolate to any $\mathbf{q}$ 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. In fact, the energy bands have already been made and preserved in files bnds.blue (LDA), bnds.red (1-shot GW) bnds.green (QSGW).

We can readily translate the bands into picture using a combination of the plbnds utility and fplot utility. Of course there are many ways to skin this cat and you may have a preferred method.

The following render the three energy bands files into the standard Questaal format for 2D arrays, which fplot can read.

echo -13,10,5,10 | plbnds -fplot -ef=0 -scl=13.6 -dat=green green
echo -13,10,5,10 | plbnds -fplot -ef=0 -scl=13.6 -dat=red red
echo -13,10,5,10 | plbnds -fplot -ef=0 -scl=13.6 -dat=blue blue


plbnds also makes a script plot.plbnds which tells fplot how to build the figure.

We’d like to plot all three bands on the same figure. To do this, modify plot.plbnds as follows:

mv plot.plbnds plot.plbnds~
echo "% char0 colg=1,bold=4,clip=1,col=.3,1,.2" >>plot.plbnds
echo "% char0 colr=3,bold=4,clip=1,col=1,.2,.3" >>plot.plbnds
echo "% char0 colb=2,bold=2,clip=1,col=.2,.3,1" >>plot.plbnds
awk '{if ($1 == "-colsy") {sub("-qr","-lt {colg} -qr");sub("blue","green");print;sub("green","red");sub("colg","colr");print;sub("red","blue");sub("colr","colb");print} else {print}}' plot.plbnds~ >> plot.plbnds  These steps modify plot.plbnds so that it will plot all three families of bands. Now do: fplot -f plot.plbnds open fplot.ps  Note: the existing lmf basis set is a little undersized for GW calculations in open systems. The new Jigsaw Puzzle Orbital basis should surmount this problem, but for the present, the best solution is to add empty sites to fill the interstices. See the Exercises for more details. ### 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 $H_0$. You can use any available $H_0$, e.g. LDA, Here we use the QSGW $H_0$ generated in the previous section. To make a good dielectric function it is necessary to modify GWinput a little. 1. The 4×4×4 k mesh used for QSGW is reasonably sufficient for k integrated properties, 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  2. 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)  3. The macroscopic dielectric function is $\varepsilon_M(\omega) = \left[\varepsilon^{-1}_{G=0,G^\prime=0}(q{\rightarrow}0,\omega)\right]^{-1}$. Because it takes an average over the unit cell, the detailed shape of $\varepsilon$ inside the unit cell is not so important. You can reduce the $L$ cutoffs significantly with minimal change in $\varepsilon_M$. Reducing it to 2 makes negligible difference in this case: lcutmx(atom) = l-cutoff for the product basis 2 2  4. Supply GWinput with the q points for which you want to calculate $\varepsilon_M(\omega)$. The code does not have the ability to take the $\mathbf{q}{\rightarrow}0$ limit analytically; you must choose a finite but small q. At the end of GWinput there are the following lines QforEPSIBZ off &lt;QforEPS&gt; 0d0 0d0 0.015d0 &lt;/QforEPS&gt;  It tells the GW code to make $\varepsilon_M(\omega)$ just a one point, close to the $\Gamma$ point. This corresponds approximately to $\varepsilon_{zz}(q=0,\omega)$. If you want to calculate $\varepsilon_{xx}$ also, add a line 0.015d0 0d0 0d0. Si is cubic and $\varepsilon_{xx}=\varepsilon_{yy}=\varepsilon_{zz}$, so it is not necessary here. You can either calculate the proper macroscopic dielectric function, $\varepsilon_M(\omega) = \left[\varepsilon^{-1}_{G=0,G^\prime=0}(q{\rightarrow}0,\omega)\right]^{-1}$, or one that neglects local fields, $\varepsilon_M(\omega) \approx \varepsilon_{G=0,G^\prime=0}(q{\rightarrow}0,\omega)$. 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  Columns 1-4 are the q-point and frequency, respectively. Columns 5-6 are the real and imaginary parts of $\varepsilon_M$; columns 7-8 are the real and imaginary parts of $\varepsilon^-1_M$. 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 $\mathrm{Im}\varepsilon^-1_M$ is zero until the first direct gap transition at 0.2516D Ry, or 3.42 eV. You can read the static dielectric constant, $\varepsilon_\infty=\varepsilon^{-1}_M(\omega{=}0)$ from EPS0001.dat. The calculated value is 9.3, which is 79% of the actual $\varepsilon_\infty$ (11.7 at 0K). This error largely originates from the neglect of ladder diagrams. Use your favorite graphics package to draw figures. To draw Re ε as a solid line and Im ε as a dotted one, do the following: $ mcx -f1p8e20.12 -r:s=1 EPS0001.dat -coll 4,5,6 -p -coll 1 -s13.6 -tog -coll 2:nc -ccat -inc 'x1>0&x1<10' > dat.lfc  # converts to standard Questaal 2D array format
$fplot dat.lfc -coll 3 -lt 3,bold=4 dat.lfc$ open fplot.ps


Compare to the neglect of local fields (superpose, use blue for no local fields)

$mcx -f1p8e20.12 -r:s=1 EPS0001.nlfc.dat -coll 4,5,6 -p -coll 1 -s13.6 -tog -coll 2:nc -ccat -inc 'x1>0&x1<10' > dat.nolfc$ fplot dat.lfc -coll 3 -lt 3,bold=4 dat.lfc -lt 1,col=0,0,1 dat.nolfc -coll 3 -lt 3,bold=4,col=0,0,1 dat.nolfc


You should see that the omission of local fields has a modest effect.

#### 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, $\varepsilon_0$. Taking the real part of $\sigma$ and scaling to appropriate units,

Using

we can write

Then

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) Add empty spheres or floating orbitals

Note: the existing lmf basis set is a little undersized for GW calculations in open systems. The new Jigsaw Puzzle Orbital basis should surmount this problem, but until it is operational there are three ways to enhance the basis

1. Add plane waves, as shown in this tutorial. You can perform 1-shot GW calculations with APWs, but not QSGW, because there is no good way to interpolate the self-energy to arbitrary k point. So it’s not a viable option here.
2. Add empty spheres — artificial sites in the interstitial with Z=0.
3. Add floating orbitals. These are like empty spheres, but have no augmentation radius. These add extra smooth Hankels to the basis centered at some interstitial point.

Until JPOs are available, the best is to use floating orbitals. Finding the right spot to place them is an art form, since there is no correct or unique point. the main guide is to identify points with large voids. The easiest way is to use blm’s automatic empty sphere finder. Do this:

$blm init.si --gmax=5 --nk=4 --nit=20 --gw --findes~float  and observe that the site file has two extra sites labelled E. --findes usually works very well, though it can add too many of them when the voids are large and symmetry is low, and you may need to manually intervene. There are several switches that can help you make better choices in such cases; see here for a modifiers to this switch. Also compare actrl.si to ctrl.si used in the tutorial. Part of the change comes from the omission of the --express switch, but also a new species was added: ATOM=E Z= 0 R= 2.221355*{(les==11)?0:1} LMX=2 LMXA=4  This line accommodates both empty spheres and floating orbitals. If variable les is 11, R will be 0. R=0 tells lmf to treat the site as a floating orbital — like an empty sphere but with no augmentation radius. les=11 is the default, since its value is set at the top of the ctrl file. It’s preferable to use floating orbitals, but there is a small catch. lmfa doesn’t now how to define envelope functions without some knowledge of the augmentation spheres, and it won’t add basis set information to the basp file. As a first step, treat E as an empty sphere, merely for the purpose of enhancing the basp file. $ cp actrl.si ctrl.si
$lmfa -vles=1 si && cp basp0.si basp.si  Note that the basp file has an extra line for species E. Repeat the tutorial with this set of (ctrl,site,basp) files, and see how much the bands change. 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


### Footnotes

1](#footnotes-and-references)) In some circumstances, e.g. to break time reversal symmetry inherent in the LDA, you need to start with LDA+U. CoO is an instance of this.