Dielectric function and other optical properties with a band code
This tutorial shows how to calculate optical and related properties for Si using either the DFT code lmf or the ASA version, lm. It has the ability to compute Im ε(q=0) in the random phase, or independent particle approximation, and also neglecting local fields. This makes execution very efficient. The GW code has much more flexibility, but it is more costly and it not covered in this tutorial. There is a related tutorial, looking at Fe. It covers similar ground also using the band code lmf, instead of the full GW machinery. It shows some extra features and uses a QSGW potential instead of a DFT one.
Preliminaries
It assumes you have completed an introductory tutorial for Si, either the full potential tutorial or steps 1 and 2 in the ASA tutorial. Your working directory should an operational input file ctrl.si, and other files generated to make the self-consistent density.
It also assumes the Questaal executables are installed somewhere in your path. This tutorial uses the ghostview utility gs view postscript files; you can substitute your favorite one.
Tutorial
Introduction
The full-potential (FP) and the Atomic Spheres Approximation (ASA) implementations of the code are executed through lmf and lm respectively, and have the capacity to perform a number of equilibrium and non-equilibrium optical and electronic calculations. This tutorial will only focus on the equilibrium calculation for optical properties and the joint density of states (JDOS); non-equilibrium modes are covered in this tutorial.
Setup for OPTICS
lmf and lm can calculate the imaginary part of the dielectric function, and some other properties. To run the optics branch, either code needs an OPTICS category and a MODE token. To see the available options, do:
lmf --input | grep -A 15 OPTICS_MODE
You should see the following:
OPTICS_MODE opt i4 1, 1 default = 0
0: do nothing;
1: generate linear eps_2
2: generate linear eps_2, spin 2
Add 20 for Second Harmonic Generation
8: eps_2 for nonequilibrium absorption. KT, IMREF are required inputs.
9: eps_2 for photoluminescence. KT, IMREF are required inputs.
-1: generate joint DOS
-2: generate joint DOS, spin 2 (LTET=3)
-3: generate up-down joint DOS (LTET=3)
-4: generate down-up joint DOS (LTET=3)
-5: generate spin-up single DOS (LTET=3)
-6: generate spin-dn single DOS (LTET=3)
-8: Joint DOS for nonequilibrium absorption. KT, IMREF are required inputs.
-9: Joint DOS for photoluminescence. KT, IMREF are required inputs.
So that we can reuse the input file in different contexts, create a reference file which will be left untouched. In the ASA context, the extra tag EFMAX is necessary so lm doesn’t take a shortcut and calculate only some of the unoccupied states. Do the following:
cat ctrl.si | sed 's/\(BZ .*\)/\1 EFMAX=99/' > ctrl.bk
In the FP context, ctrl.si as given will assign BZ_NKABC to 4. We will want to vary this number and to do this use the stream editor sed to make the change
echo '% const nk=4' > ctrl.bk
sed 's/nkabc= *4/nkabc={nk}/' ctrl.si >> ctrl.bk
Whether ASA or FP, make sure you have followed the steps to make the density self-consistent.
Dielectric function
The imaginary part of the dielectric function, Im ε, is computed with OPTICS_MODE=1. Create the input file as follows:
echo OPTICS MODE=1 NPTS=1001 WINDOW=0 1 LTET=3 > ctrl.si
cat ctrl.bk >> ctrl.si
When either lmf or lm is executed, it will generate Im ε. Both setups used a fine enough k mesh for charge self-consistency, but for optics a much finer mesh is needed to get a clean figure. Run lm or lmf with extra k points:
lm si --quit=band -vnk=16 or lmf si --quit=band -vnk=16
This creates and stores the longitudinal parts of Im ε in file opt.si. (The longitudinal elements of ε are εxx, εyy, εzz i.e. the diagonal part of the tensor εij). Si has cubic symmetry, so all three columns should be the same.
Inspect Im ε visually. Note that the first column is the frequency, in Ry. The following draws creates a figure for Im ε with the abscissa in eV
fplot -ab 'x1*13.6' opt.si
gs fplot.ps
There should be a sharp rise near 2.7 eV. This is the DFT direct gap in Si (experimentally it is closer to 3.3 eV). (There is an indirect gap at 1.2 eV, but those transitions require a phonon assist which this code does not have the ability to include).
Re ε can be obtained from Im ε using a Kramers-Kronig transformation. You can find matlab or mathematica codes on the web to do this, but they are not part of Questaal.
The Questaal package has a very primitive facility to make Re ε from Im ε, using an old “dusty deck” program, called kram. It requires an input a file, named kram.man, along with Im ε. Cut and paste the following box into kram.man:
in
999 !NMAX this is number of energy points
1 0.1 !ISMO,GA(EV) 1 means yes: smooth, GA is gaussian width
0 0.607 469. 200. !IDRUDE=1 mean drude inclued,V2N,VOL,RO (drude params)
5. !Theta (degrees), from surface (so 5 degree is near glancing)only meaningfull if NQQQ=3 see below
outfile
2 !NQQQ(1:E2,2:E1,3:REF,4:EELS,5:E2(i),6:SIG(i)
! this selects what to do
1 is epsilon2 (menaning just convert E to eV and possibly smooth
2: epsilon1
3: reflectivity
4: EELS gives Im (-1/eps)
5: not sure what this one gives!
6: conductivity sigma eps=1 + 4 pi i/omega sigma(omega)
The first line tells kram to reads from file in. We must strip the first line from opt.si, as kram doesn’t recognize it. We will do this and redirect the output to file in. Then we can run kram directly:
grep -v cols opt.si > in
kram
It dumps Re ε into outfile (the name is specified in line 6 of kram.man). (Note that kram converts the abscissa from Ry to eV, and actually returns Re ε−1). Visually inspect outfile:
fplot opt.si
gs fplot.ps
You should find that the static ε (called ε∞, because the true static dielectric constant includes the contributions from nuclear displacements) is about 10. This is a little smaller than the experimental constant of 11.7. It is understimated because the RPA leaves out electron-hole attraction in the polarizability. Note also that Re ε has a zero at around 4 eV. Re ε crossing zero corresponds to a plasmonic excitation.
Joint Density-of-States
The joint Density-of-States (DOS) is very similar to Im ε, except that the optical matrix element is left out. It also has an extra ability, to join unoccupied states at a wave number shifted from the occupied ones. The normal j-DOS (vertical transitions) proceeds exactly as for Im ε, except for a different OPTICS_MODE. For the ASA do
echo OPTICS MODE=-1 NPTS=1001 WINDOW=0 1 LTET=3 > ctrl.si
cat ctrl.bk >> ctrl.si
lm si --quit=band -vnk=16 or lmf si --quit=band -vnk=16
Data is written to jdos.si. Draw a figure to confirm a sharp onset at 2.7 eV.
fplot -ab 'x1*13.6' jdos.si
gs fplot.ps
Joint Density-of-States, transitions with offset q
Si has an indirect gap near q=(0,0,0.85) in units of 2π/a, so it would be of interest to see the joint DOS for transitions with this q. (The proper dielectric response at q=(0,0,0.0) will include these transition to this q, cannot be calculated properly because the dielectric function it involves an electron-phonon matrix element not built into Questaal.)
Transitions can be calculated for any q on the discrete mesh, so we must find the mesh point closest to q=(0,0,0.85). We can express q as a linear combination of reciprocal lattice vectors, Qlat = (-1,1,1, 1,-1,1, 1,1,-1), as follows:
mcx '-array[3,3]' -1,1,1,1,-1,1,1,1,-1 -i '-array[3,1]' 0,0,0.85 -x -s16
Expressed as multiples of Qlat, q=(0.425,0.425,0). On a discrete mesh of 16 divisions, this translates into mesh point (6.8,6.8,0) as the preceding mcx calculation will show. The mesh point of closest approach is then (7,7,0). To calculate the joint DOS for this q, add a line IQ=7,7,0 to the OPTICS category, e.g.
echo OPTICS MODE=-1 NPTS=1001 WINDOW=0 1 LTET=3 > ctrl.si
echo ' IQ=7,7,0' >> ctrl.si
cat ctrl.bk >> ctrl.si
lm si --quit=band -vnk=16 or lmf si --quit=band -vnk=16
Notice that the output has this block:
OPTINQ: JDOS (++) using occ bands (1,26) and unocc (1,26)
1001 points in energy window (0,1) Efermi=0.185517
q for unocc states: iq = 7 7 0, q = 0.000000 0.000000 0.875000
Reduce bands to range (1,14)
It confirms that iq = 7 7 0 corresponds to a q transfer of 0.875, as we required.
Draw the figure and you should now see the onset at about 0.5 eV, corresponding to the indirect gap of Si.
fplot -ab 'x1*13.6' jdos.si
gs fplot.ps
If you want ε(ω,q) for any q, you must use the optics branch of the GW code.
Limit occupied and unoccupied bands
It is possible to perform any of the optics mode calculations described above for a restricted number of bands. This can greatly speed up the calculation and allow for isolation and identification of individual band contributions. To restrict the bands involved in the calculation simply provide a range of values for occupied and unoccupied bands through OPTICS_FILBND and OPTICS_EMPBND respectively. Below is an example of an optics category which calculates the contribution to Im from the highest three valence bands and the lowest two conduction bands for Si:
echo OPTICS MODE=1 NPTS=1001 WINDOW=0 1 LTET=3 > ctrl.si
echo ' FILBND=3,4 EMPBND=5,6' >> ctrl.si
cat ctrl.bk >> ctrl.si
lm si --quit=band -vnk=16 or lmf si --quit=band -vnk=16
It is written to file opt.si, as before. This data is the same as the original optics calculation, but decays a little more rapidly.
Resolve into contributions from individual pairs
Because optics are approximated by the independent-particle approximation, Im ε can be resolved by contributions from individual (occ,unocc) pairs. To turn on this feature, simply add PART=1 to the OPTICS category, e.g.
echo OPTICS MODE=1 NPTS=1001 WINDOW=0 1 LTET=3 > ctrl.si
echo ' FILBND=3,4 EMPBND=5,6' >> ctrl.si
echo ' PART=1' >> ctrl.si
cat ctrl.bk >> ctrl.si
lm si --quit=band -vnk=16 or lmf si --quit=band -vnk=16
Now the result is written into file popt.si. This file has 12 columns, in addition to the first column (which is the energy). Columns 2-5 correpond respectively to transitions between (3,5), (4,5), (3,6), (4,6) pairs, for Im εxx; columns 6-9 are for Im εyy and columns 10-13 are for Im εzz.
You can confirm that four resolved Im εxx sum to the total:
mc popt.si -e2 x1 x2+x3+x4+x5 opt.si -coll 1,2 --
Resolve into contributions by wave number
Im ε can be also resolved by wave number k. To accomplish this, use OPTICS_PART=2:
echo OPTICS MODE=1 NPTS=1001 WINDOW=0 1 LTET=3 > ctrl.si
echo ' FILBND=3,4 EMPBND=5,6' >> ctrl.si
echo ' PART=2' >> ctrl.si
cat ctrl.bk >> ctrl.si
lm si --quit=band -vnk=16 or lmf si --quit=band -vnk=16
Now there are 135×3 channels, since there are 135 irreducible k points, so popt.si gets pretty unwieldy. You may wish to write the results to a binary file. There is also an optics editor to help you analyze results; see this page.
Resolve into contributions by pairs and by wave number
Finally, you can resolve Im ε by both pairs and wave number. Follow the same procedure, but with OPTICS_PART=3.
Further ways to resolve DOS or Im ε
You can also project the eigenfunctions onto particular channels, which gives you a partial DOS from those channels. You can further combine this with the PART options above. For a tutorial and demonstration, see this tutorial.