# 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 QS*GW* 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

*ε*and columns 10-13 are for Im

_{yy}*ε*.

_{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.