# Equation of states tutorial {#sec:eos}

## Introduction

In this tutorial, you will see how to compute the lattice constant of
diamond using total energy calculations. We will use a Birch--Murnaghan
equation of state fit to extract the equilibrium volume. The
Birch--Murnaghan equation of state reads

$$E(V) = E_0 + \frac{9V_0B_0}{16}
\left\{
\left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^3B_0^\prime +
\left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^2
\left[6-4\left(\frac{V_0}{V}\right)^\frac{2}{3}\right]\right\}.$$

## Write an input file template

We will perform 7 total energy calculations around the experimental
value of 3.57 A. 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'
aa = 3.57 * xfacx / 100
```

There is a tag `xfacx` that will be replaced by numbers, and hence we
call this a template. For example, setting the tag to `100` will make
the lattice constant `3.57 * 100 / 100 = 3.57`. The first two lines set
the name of the result file and standard output log file respectively.
The third line sets the lattice constant. As for the rest of the file,
you should define the structure and set the resolution and k-sampling to
sensible values. If the accuracy parameters are not stringent enough,
the fit will not be of good quality.

## Perform the calculations

As mentioned above, we will perform 7 total energy calculations around
the experimental value of 3.57 A. To generate the 7 corresponding input
files, one may run the script `genall.sh <genall.sh>`{.interpreted-text
role="download"}. The important thing here is `{94..106..2}` which means
that the tag will be replace by `94, 96, 98, 100, 102, 104, 106`.
Depending on the material, this range may be too large or too narrow.
You should adjust it so that the energy varies by at least a few meV.
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
simultaenous 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_eos.m <read_eos.m>`{.interpreted-text role="download"}. It will
list the files that match the pattern `results/scf_f*.mat`, there should
be 7 in total. Then it will call the function `eosfit` to do the fit.
Finally, it will print

``` shell
v_eq = 5.647781 ang^3/atom
```

From that we can deduce the lattice constant for the conventional
diamond lattice. Multiply by 8 and take the cubic root, and get `3.56`
A. You may also plot the EOS as follows

``` matlab
figure(); hold on
plot(avol,anrg)
xlabel('volume (A^3/atom)')
ylabel('energy (eV/atom)')
saveas(gcf,'eos_diamond.png')
```

which gives

<figure>
<img src="./eos_diamond.png" alt="./eos_diamond.png" />
<figcaption>Diamond equation of states</figcaption>
</figure>
