# Valence band maximum tutorial {#vbm_tutorial}

## Introduction

In this tutorial, you will see how to compute the valence band maximum
of diamond. This VBM is obtained invoking Janak's theorem

$\frac{d E(n_i)}{dn_i} = e_i$

where $n_i$ is the occupation number of the highest occupied state $i$
with the eigenvalue $e_i$. This is calculated by finite difference, i.e.
by removing a small fractional charge from the diamond unit cell.

$e_i = \frac{E[Q=N] - E[Q=N-\delta N]}{\delta N}$

where $\delta N = 10^{-3}$ is usually a good value.

## Write an input file template

We will perform 4 total energy calculations to obtain 3 VBM values. We
first write an input file template like
`scf_base.input <scf_base.input>`{.interpreted-text role="download"}.
Note the lines

``` shell
info.savepath = './results/scf_fxfacx'
info.outfile = 'fscf_fxfacx.out'
atom.nvalence = 2 * 4 - 1 / xfacx
```

The tag `xfacx` will be replaced by the inverse fraction of removed
electrons. For example, setting the tag to `100` will set the valence to
`2 * 4 - 1 / 100 = 7.99`, removing 0.01 electron. The first two lines
set the name of the result file and standard output log file
respectively.

## Perform the calculations

As mentioned above, we will perform 4 total energy calculations. To
generate the 4 corresponding input files, one may run the script
`genall.sh <genall.sh>`{.interpreted-text role="download"}. The tag
values are set to `10, 100, 1000, 10000`. Our template uses the
pseudopotential file `C_SZDP.mat <../C_SZDP.mat>`{.interpreted-text
role="download"} which should be at the root of the `tutorials`
directory. One can perform the calculations one by one, or run them in
parallel using the script `runall.sh <runall.sh>`{.interpreted-text
role="download"}. Before using `runall.sh`, modify the script
considering that `N` is the number of simultaneous calculations and
`NCORE` is the total number of cores you want to use.

## Determine the equilibrium volume

The last step is to extract the equilibrium volume from the equation of
states. This is done with the script
`read_vbm.m <read_vbm.m>`{.interpreted-text role="download"}. It will
list the files that match the pattern `results/scf_f*.mat`, there should
be 4 in total. Then it will compute the VBM values using finite
differences. Finally, it will print

``` shell
vbm = 10.126067 eV
```

The smallest difference value is printed, but make sure that it is
converged. This is readily seen on a plot of the VBM as a function of
difference.

``` matlab
figure()
semilogx(x,vbm); hold on
xlabel('\Delta N')
ylabel('\Delta E / \Delta N (eV)')
saveas(gcf,'eos_vbm.png')
```

which gives

<figure>
<img src="./eos_vbm.png" alt="./eos_vbm.png" />
<figcaption>VBM versus occupancy difference</figcaption>
</figure>
