Silicon – full frequency

In this session we will calculate the quasiparticle band structure of silicon
using the LDA and GW approximations, with full-frequency integration.

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 GW calculation on silicon with full-frequency integration.
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.

Instructions
————

Mean-field calculation

As we did in the Silicon-generalized-plasmon-pole-model example, we first perform all the mean-field calculation required 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. Follow these steps:
    • 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
      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 2.2-wfnq 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 “WFNq“ by the Epsilon code. Follow these steps:
    • 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
      output “kgrid.out“ and the more verbose log “kgrid.log“.
      – 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.
    • 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 4-bandstructure directory. In this directory we perform a non-self-consistent (NSCF) mean-field calculation to generate all of the Bloch states needed for the interpolation of GW calculations on to a k-path. The resulting wavefunction file will then be read as “WFN_co“ by the Inteqp code. Follow these steps:
    • 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-full-frequency

Now we are ready to perform GW calculations. Create and go to the “2-bgw/“ directory, where we 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. Follow these steps:
    • 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 1/epsilon^{-1}_{0,0}(0) and epsilon_{0,0}(0)?
      Why is there a difference? (Hint: find these values by searching for
      “Head“.) Make a plot of 1/epsilon^{-1}_{0,0}(q) and epsilon_{0,0}(q)
      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.
      Try to understand their meaning and whether there is any cause for concern.
    • Look at “EpsDyn“. Plot it using your favorite tool (e.g. gnuplot). Does the structure of the file make sense?
  2. Create and go to 2-sigma directory. In this directory we calculate the GW quasiparticle energies for silicon using the Sigma code from BerkeleyGW. Follow these steps:
    • 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 the file “spectrum.dat“ – this shows you an energy-resolved spectrum from the full-frequency calculations. Look at the imaginary parts of the self-energy. Does the structure make sense?
    • Plot the real and imaginary parts of the self-energy for a couple of bands and k-points. What can you learn?
  3. Create and go to 2-sigma/bandstructure. In this directory we calculate an interpolated quasiparticle bandstructure for silicon using the Inteqp code from BerkeleyGW. Follow these steps:
    • 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 data

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