# Dielectric constant tutorial {#dielectric_tutorial}

## Introduction

In this tutorial, you will see how to compute the dielectric constant of
diamond using an external electric field. In essence, the electrons
react to the introduction of an external field by screening it. By
looking at the variation of the effective (Kohn-Sham) potential we can
determine by how much. The ratio of the external potential to the
effective potential variation is the dielectric constant.

## The calculations

We introduce a sawtooth potential, which has the virtue of being
periodic and linear away from the bends, making the analysis
straightforward.

We first perform the calculation defined in
`scf_bulk.input <scf_bulk.input>`{.interpreted-text role="download"},
which will yield the reference effective potential. Note the lines

``` shell
option.replicateCell = [1,1,8]
option.savePotential = 1
```

The input file `scf_bulk.input` defines a conventional cubic unit cell
of diamond. But the keyword `option.replicateCell` tells RESCU to set up
a $1 \times 1 \times 8$ supercell calculation. Then
`option.savePotential` is set to `1` so that the potential terms are
saved before exiting the program.

We then perform the calculation defined in
`scf_vext.input <scf_vext.input>`{.interpreted-text role="download"},
which will yield the screened effective potential. Note the lines

``` shell
potential.external = @(u,v,w) sawtooth(0.5,u,v,w,3)
rho.in = 'results/scf_bulk'
```

The keyword `potential.external` allows introducing an external
potential as a function handle. Here sawtooth is a user defined function
whose arguments are: amplitude, u-grid, v-grid, w-grid and dimension. It
should be found in `src`. The amplitude is in units of eV. The u-grid,
v-grid, w-grid are fractional coordinate grids. The dimension is `3`,
meaning that the sawtooth function runs along the third lattice vector
direction. Finally, we set `rho.in` to `results/scf_bulk` to restart
from the bulk results.

## Determine the dielectric constant

The last step is to extract the dielectric constant from the potentials.
This is done with the script `read_eps.m <read_eps.m>`{.interpreted-text
role="download"}. Basically, the script does the following

1.  Load the potentials.
2.  Perform a linear fit of the external potential in the center of the
    linear potential region.
3.  Perform a linear fit of the effective potential in the center of the
    linear potential region.
4.  Determine the dielectric constant as the ratio of the two.

Finally, it will print

``` shell
dielectric constant = 5.410509
```

You may also plot the potentials and fit as follows

``` matlab
figure(); hold on
z = w * avec(3,3);
tmp = squeeze(mean(mean(reshape(vext, fgn),1),2)); 
plot(z,tmp, '--', 'DisplayName', 'v_{ext}'); hold on
tmp = squeeze(mean(mean(reshape(vel - velb, fgn),1),2)); 
plot(z,tmp, '--', 'DisplayName', '\Delta V'); hold on
vext = sawtooth(0.5/p1(1)*p2(1),u3,v3,w3,3);
tmp = squeeze(mean(mean(reshape(vext, fgn),1),2)); 
plot(z,tmp, '--', 'DisplayName', 'v_{fit}'); hold on
xlabel('z (bohr)')
ylabel('potential (Ha)')
legend(); grid on
saveas(gcf,'vext_diamond.png')
```

which gives

<figure>
<img src="./vext_diamond.png" alt="./vext_diamond.png" />
<figcaption>Potential screening in diamond</figcaption>
</figure>
