Every thing we’ve done so far, has in fact been for systems at effectively zero temperature. And even at zero temperature, we haven’t taken any account of the zero-point motion of the atoms and have treated them as purely classical objects residing at fixed positions.

In actuality, there will be an effect from the zero-point motion of the atoms, and as temperature increases the atoms will vibrate with larger amplitude and this can lead to changes in many properties of a material with temperature.

In this lab we’ll be looking at a few ways of calculating this.

This is quite a challenging problem, and is the focus of active current research. We’ll be looking at how one might approach this for some simpler examples, but there are many possible ways to calculate these kinds of properties.

Our approach will be to use the type of DFT and DFPT calculations you’ve already seen how to do, and spend more time analysing the output to produce materials properties.

This will involve some reasonably serious numerical calculations, of the kind that would be quite difficult to do without using some form of mathematical software. In this lab, you’ll be using jupyter ipython notebooks. These allow us to mix text/latex descriptions with python code.

We’ll be mainly relying on the following libraries:

- numpy
- This allows us to easily work with arrays of data.

- scipy
- This gives us many numerical analysis tools.

- matplotlib
- This allows us to easily generate plots.

I expect many of you may not be familiar with python, but if you’ve used Matlab or Mathematica, you should find the process here somewhat similar, but with slightly different commands.

Jupyter and Python are freely licensed, so if you want to install them on your own computer, you can either follow the above links or you can get a complete package from anaconda which will give you everything you need on Windows, Mac or Linux. There is extensive documentation and usage examples at the above links.

To begin with, make a copy of the directory for this lab as usual, and in your
copy type `jupyter notebook`

. This will output some text to your terminal,
and you may see some error messages associated with running this over x2go,
but in a minute or so it should open a window on your browser with a graphical
interface to the notebook. Each of the folders contains a notebook,
suffixed `.ipynb`

associated with analysis of the DFT calculations in that
folder.

To evaluate the code in a jupyter cell press `ctrl+enter`

, or press
`shift+enter`

to evaluate the cell and automatically move to the next cell.

We’ll be using many numerical functions, and if you want information on what a
function call we’re using is doing, you can enter `help(functionname)`

in a
cell and it will output a detailed description of the input, output and method
used once you evaluate the cell.

jupyter will take over the terminal you start it in, so it’s best to open a second terminal to do your calculations in. If you like, you can start a second tab in the same terminal application (if you’re using the default) by pressing

`ctrl+shift+t`

.

## Average Bond Length vs Vibration Amplitude of a Molecule

To give you an initial idea of how this works, let’s look our H2 molecule
again. This is in the `01_H2_bl`

directory.

Here we have set up a template file, `H2_base.in`

, that we’ll use to calculate
the total energy of a H2 molecule as a function of the separation between the
two atoms. The script `EvBL.sh`

will perform the calculation, incrementing
the space between the atoms each time and generating a new input and output
file, accumulating the energy vs bond length in `etot_v_bl.dat`

.

Run this now, and plot it in gnuplot to ensure it looks reasonable.

If it all looks OK, let’s go back to the jupyter notebook you opened earlier
in your browser and we’ll continue the analysis there. Click on the
`01_H2_bl`

directory, and then click on
`BLvT.ipynb`

to open the notebook. You’ll need to
execute the various cells containing code in the order they appear (click on
them to highlight them, and press `ctrl+enter`

) to perform the
analysis.

## Thermal Expansion of a Solid

Next we’ll look at the thermal expansion of a solid. We’ll do this in three parts:

- Examine the behaviour of the DFT total energy as a function of volume.
- Examine the effect of including the vibrational energy to calculate the free energy at zero temperature and find the effect on the equilibrium volume,
- Calculate the volume that minimizes the free energy as a function of temperature. From this we can obtain the thermal expansion.

### Step 1 - DFT Total Energy vs Volume

To perform step 1, take a look in the `02_Si_V`

folder. This is
fairly similar to the H2 calculation we performed earlier and involves
calculating the total energy of silicon as a function of volume. The
calculation is set up using a similar scripted loop as you’ve seen before.

Run this now, and plot the resulting `etot_v_vol.dat`

in gnuplot to
ensure it the calculation worked correctly.

If it all looks OK, open the jupyter notebook in browser and we’ll continue the analysis there.

### Step 2 and Step 3 - Free Energy vs Volume and Temperature

We’ll be using the same DFT output for steps 2 and 3, which will be (in
fact this has already been calculated and made available to you)
generated in the `03_Si_F`

folder. The key extra thing we
need for these steps is the phonon density of states at several different
volumes above and below the equilibrium. There is a script set up in this
directory that was used to produce the phonon DOS at 6 different volumes.

Take a look at the script `FvV.sh`

and ensure you understand it. It’s an
expanded version of the type of script you’ve seen before. Recall, we needed
to combine several different codes from the quantum espresso package to
calculate the phonon DOS. This script uses a template file for the input for
each of those codes, and along with modifying the cell volume in the
SCF calculation also sets a prefix for the calculation index so that the
various output files are named differently for each volume. Either
1 or 2 values will be modified from the template in each input file.

If you recall, it took around 5 minutes to calculate the phonon DOS for a
single volume last time. Doing this 6 times would take half an hour. For this
reason, I’ve included the dos files that would be generated. **You don’t
need to run the calculation here and can proceed to the analysis in the
jupyter notebook in this directory.**

## Summary

In this lab you have seen:

- For a molecule we used DFT calculations of the energy vs bond length to show how the average separation increases with amplitude of oscillation.
- For a solid we used
- DFT calculations of the total energy vs volume with a fit of an equation of state.
- DFPT calculations of the phonon mode density of states vs volume with a fit of the vibrational energy vs volume.
- This allowed us to get both zero temperature and finite temperature corrections to the equilibrium volume, coming from the vibrational energy and entropy.