Silicon – generalized plasmon pole model

In this session we will calculate the quasiparticle band structure of silicon
using the LDA and GW approximations using a plasma-pole model. Notice that we use PARATEC as our choice of mean-field code here but BerkeleyGW has also been interfaced with Quantum Espresso, ABINIT, Octupus, etc.

Goals
—–

The basic goals are the following:

1. Understand the basic workflow of BerkeleyGW, and the relation between the k-grids, wavefunctions, Epsilon, Sigma, and Inteqp.
2. Run a basic GW calculation on silicon.
3. Construct an interpolated bandstructure via scissors parameters and Inteqp.

The stretch goals are:

1. Compare your Sigma GW results with Hartree-Fock and/or static COHSEX.
   What inputs are no longer necessary? How do the results compare?
2. Modify the example for GaAs and repeat each step of the calculation.
3. Adapt the calculation to use half-shifted k-grids.
4. Adapt the calculation to use Quantum ESPRESSO or EPM as a mean field.

Instructions
————

Mean-field calculation

In the first step we perform all the mean-field calculation needed in the GW calculations. Create and go to the “1-mf/“ directory, where we perform all the mean-field calculations.

1. Create and go to 1-scf directory, In this directory we perform a standard SCF calculation and obtain a charge density which will be used to generate all the other WFNs. We use a standard Monkhorst-Pack k-grid that converges the charge density, and no input flags pertaining to BerkeleyGW need to be used for PARATEC here.
   • Start the calculation with “./01-calculate_scf.qsub“
   • Take a look at “input“ for PARATEC. What functional, cutoff, and k-grid are being used? Look for key input flags described in the manual section MeanField/PARATEC/README (most are not present at this stage).
   • Look briefly at the “OUT“ file when the run finishes. Did the run succeed? Find the list of k-points in the “OUT“ file and “SCF_KPOINTS“.
2. Create and go to 2.1-wfn directory, In this directory we perform a non-self-consistent (NSCF) mean-field calculation to generate all of the Bloch states needed for the subsequent GW calculations. The resulting wavefunction file will then be read as “WFN“ by the Epsilon code, and “WFN_inner“ by the Sigma code.
   • Run “./02-calculate_wfn.qsub“ to calculate the wavefunctions with PARATEC.
   • Take a look at the input file (“kgrid.inp“) for “kgrid.x“, the executable that generates our k-grids for NSCF calculations. Compare this input file with the documentation available in the code. What is the meaning of each field? How many k-points do you expect to be generated? Find the FFT grid in the PARATEC output file in “1-scf“ and check that it is correct in “kgrid.inp“. Why do we need this? What happens if you use a different value (e.g. an odd one)?
     How does the kgrid used here compare to that in “1-scf“?
   • Run “./01-get_kgrid.sh“ to generate the k-grid. Take a look at the
     output “kgrid.out“ and the more verbose log “kgrid.log“.
     – What is the space group?
     – What are the symmetry elements of the crystal?
     – What k-grid was generated?
     – How many k-points are there in the full and symmetry-folded Brillouin zone?
     – See what symmetry operation was used to relate two points in the full BZ and reduce them to only one in the final result.
   • Take a look at the “input“ file for PARATEC. How does it compare to that
     of “1-scf“? Find the key input flags described in the manual in
     MeanField/PARATEC/README.
     What does an NSCF run mean, and why do we do it instead of just another SCF? Find the “k_grid“ parameters and compare to those in “kgrid.inp“. PARATEC is able to perform the same construction as in “kgrid.x“, but for other codes we would need to copy “kgrid.out“ to the input file. Should you be writing a real or complex wavefunction? Why?
   • When the job is done, compare the k-points generated by PARATEC (listed in the “OUT“ file) to those generated by “kgrid.x“. Why do the k-points have the weights they do?
   • Compare the symmetry elements reported by PARATEC to those in
     “kgrid.log“. Run “wfn_rho_vxc_info.x WFN“ and take a look at what information is provided. Compare the symmetry elements here to those from PARATEC (from which they are derived). Relate the parameters of the PARATEC calculation to what is in the WFN file.
   • Inspect the “RHO“ and “VXC“ files with “wfn_rho_vxc_info.x“. What
     information is here compared to for “WFN“? Look at the ASCII file
     “vxc.dat“. How do the matrix elements here relate to “WFN“ and “VXC“? Are they real or complex? Why?
3. Create and go to directory 2.2-wfnq. In this directory we perform a non-self-
   consistent (NSCF) mean-field calculation to generate all of the Bloch states
   needed for the subsequent GW calculations. The resulting wavefunction file will
   then be read as “WFNq“ by the Epsilon code.
   • What k-grid was generated?
   • How many k-points are there in the full and symmetry-folded Brillouin
     zone?
   • Try to match up some points between this “kgrid.out“ and that of “2.1-wfn“: each point in “WFNq“ corresponds (with perhaps applying a symmetry) to one in “WFN“, plus the extra shift.
   • Run “./02-calculate_wfnq.qsub“ to calculate the wavefunctions with PARATEC.
   • Take a look at the input file (“kgrid.inp“) for “kgrid.x“, the
     executable that generates our kgrids for NSCF calculations. How do the
     parameters here differ from in “2.1-wfn“? Why do we need a shifted grid for this calculation?
   • Run “./01-get_kgrid.sh“ to generate the k-grid. Take a look at the
     output “kgrid.out“ and the more verbose log “kgrid.log“.
   • Take a look at the “input“ file for PARATEC and compare to the key input flags described in the manual in MeanField/PARATEC/README. How does the file compare to that of “2.1-wfn“? Do we have any unoccupied bands? Why?
   • When the job is finished, compare the k-points generated by PARATEC (listed in the “OUT“ file) to those generated by “kgrid.x“.
4. Create and go to directory 4-bandstructure. In this directory we perform a
   non-self-consistent (NSCF) mean-field calculation to generate all of the Bloch
   states needed for the subsequent interpolation of GW calculations to specified
   k-path. The resulting wavefunction file will then be read as “WFN_co“ by the
   Inteqp code.
   • Run “./01-calculate_wfn.qsub“ to calculate the wavefunctions using PARATEC.
   • Take a look at the bandlines commented out at the end of “input“, and compare to the “KPOINTS“ file generated from them. Identify what paths are being taken with respect to a diagram such as http://www.iue.tuwien.ac.at/phd/dhar/node18.html.
   • Why don’t we run “kgrid.x“ for this? Should the points be reduced by symmetry? How many bands are we using?
   • Take a look at the output file. Can you find the list of k-points?
     Use “wfn_rho_vxc_info.x“ to inspect the k-points in the “WFN“ file.
     What is listed for the k-grid, shifts, and k-weights in the file? Are they
     meaningful here?

GW calculation-GPP model

For the following part of the examples, create and go to the “2-bgw/” directory, where we start to perform GW calculations.

1. Create and go to 1-epsilon directory. In this directory we calculate the RPA dielectric matrix of silicon using the Epsilon code from BerkeleyGW.
   • Run the epsilon calculation using “./01-run_epsilon.qsub“.
     • Study the input file “epsilon.inp“ and compare to the input file
       description “Epsilon/epsilon.inp“ in the manual.
     • What files have been linked into this directory? How are they used?
   • The dielectric function epsilon(q) must be calculated for all q-points
     such that q = k – k’. For |q| /= 0, all k-points are coming from “WFN“. In
     this case, which are the possible q-points for epsilon(q)? Use “kgrid.x“ (and
     your own editing) to generate the list of q-points for “epsilon.inp“.
     Why do we need to specify q-points? What q-shift are we using?
   • Use “degeneracy_check.x“ to determine what numbers of bands
     are acceptable with respect to degeneracy from “WFN“. Does “WFNq“ matter? Is the number listed in “epsilon.inp“ an acceptable number? What happens if you change it to another value? How many bands are available in “WFN“ and “WFNq“?
   • When the job is finished, take a look at “epsilon.out“.
     What are the values of and ?
     Why is there a difference? (Hint: find these values by searching for
     “Head“.) Make a plot of and against |q|. (Hint: look for the section with “Q0 and |Q0|“.) Is it following the expected behavior as q -> 0? What do these values mean physically? Which is “the dielectric constant“ as you would measure with a capacitor?
   • Take note of the two files that were produced: “epsmat.h5“ and “eps0mat.h5“. What do they mean? Which q-points should be in each one?
   • Look at “epsilon.out“ for any warnings (marked “WARNING“) that may have occurred.
2. Create an go to 2-sigma directory. In this directory we calculate the GW quasiparticle energies for silicon using the Sigma code from BerkeleyGW.
   • Run the sigma calculation using “./01-run_sigma.qsub“.
   • Study the input file “sigma.inp“ and compare to the input file
     description “Sigma/sigma.inp“ in the manual.
   • What files are linked into this directory? How are they used?
   • Which “kgrid.x“ run should you use to populate the list of k-points in “sigma.inp“? What are these k-points used for? What happens if you add or remove some?
   • Use “degeneracy_check.x“ to determine what numbers of bands are
     acceptable with respect to degeneracy from “WFN_inner“. Is the number
     listed in “sigma.inp“ an acceptable number? How many bands are available
     in “WFN_inner“? How do the degeneracy numbers compare to those needed for Epsilon? Why?
   • When the job is finished, take a look at the output file “sigma_hp.log“.
     Compare the k-points and bands requested to what is written in the results.
     What does each column mean? Which one is the best estimate of the quasiparticle energies? What are the values of the quasiparticle renormalization factors?
   • Look at “sigma.out“ for any warnings (marked “WARNING“) that may have occurred. Try to understand their meaning and whether there is any cause for concern.
3. Create and descend into the “bandstructure“ directory to calculate an interpolated quasiparticle bandstructure with the inteqp code.
   • Run the inteqp calculation using “./01-run_inteqp.qsub“.
   • Study the input file “inteqp.inp“ and compare to the input file
     description in “BSE/absorption.inp“ (sections pertaining to inteqp). Why do we need to interpolate rather than just directly run Sigma on these k-points? How is the number of coarse points listed determined? Why do we want to use symmetries on the coarse grid but not the fine one?
   • What files are linked into this directory? How are they used?
   • Use “degeneracy_check.x“ to determine what numbers of conduction
     and valence bands are acceptable with respect to degeneracy from “WFN_co“ and “WFN_fi“. Are the numbers listed in “inteqp.inp“ acceptable? How many bands are available in “WFN_co“ and “WFN_fi“? How do the degeneracy numbers compare to those needed for Epsilon or Sigma? Why?
   • When the job is finished, take a brief look at “inteqp.out“. What warning
     do you see at the end? Look at “dcmat_norm.dat“ and “dvmat_norm.dat“ to see which values are being flagged as less accurate. How could you improve the accuracy of the interpolation for these energy levels?
   • Plot the resulting bandstructure.
     • First, load these modules:
       “module load matplotlib numpy python“
     • Then, “python plot_bandstructure.py“.
       This will produce “Si_bands.eps“. View with “display Si_bands.eps“.
   • Using “Si_bands.eps“ and “bandstructure.dat“, find the energy and location of the the direct gap and indirect gap, for LDA and GW. What is the valence band-width for LDA and GW? How do the band curvatures (effective masses) and band topology compare between LDA and GW?
   • Use “qp_shifts.py ../sigma_hp.log“ to produce “cond“ and “val“
     files containing a list of MF eigenvalues (column 1) and quasiparticle
     corrections (E_qp – E_MF, column 2) for conduction and valence bands,
     respectively, in a format suitable for plotting. Make a graph of
     quasiparticle corrections as a function of initial energy. Over what range
     of energies around the gap would a linear fit be appropriate? Why is the
     relation different for conduction and valence? Why is it non-linear farther
     away?
   • Perform a linear fit on the linear regime, separately for “cond“ and “val“.
     For example with gnuplot, to fit between energies E_1 and E_2, you would do:gnuplot > f(x)=a*x+b;a=1;b=1 # a,b initial guesses, doesn’t matter
     gnuplot > fit [E_1:E_2] f(x) ‘cond’ via a,b
     [output giving you the fit parameters]
     gnuplot > plot f(x), ‘cond’ # let you see the fitted line with dataHow good is the fit? How well does it describe the QP corrections farther from
     the gap? You will later use these fits in the calculation of the optical
     absorption, as a stretch goal in a later example.
   • Do a fit also on the whole range of values (just remove the “[E_1:E_2]” part). How good are the fits? Apply these linear relations to “bandstructure.dat“ and compare the results to the more rigorous (but time-consuming) calculation from inteqp. You can also add the inteqp “bandstructure.dat“ values to “cond“ and “val“ and see how the points there compare to the ones from Sigma on the regular grid.