Exercise 4: Lattice Expansion in the QuasiHarmonic Approximation
Estimated total CPU time: 35 min
Warning
In the following exercises, computational settings including the reciprocal space grid (tag k_grid
), the basis set, and supercell's size, have been chosen to allow for a rapid computation of the exercises in the limited time and within the CPU resources available during the tutorial session. Without loss of generality, these settings allow to demonstrate trends of the lattice dynamics of materials.
In the production calculation, all computational parameters should be converged.
In this exercise, you will:
 Perform phonon calculations in supercells with different volumes.
 Learn how to use the harmonic vibrational free energy to determine the lattice expansion.
We are going to inspect how the thermal motion of the atoms at finite temperatures can lead to an expansion (or even a contraction) of the lattice. For an ideal harmonic system, which is fully determined by the dynamical matrix \(D_{IJ(\mathbf{q})}\) defined in Eq. (3), the Hamiltonian (cf. Eq. (1)) does not depend on the volume. This also implies that the harmonic Hamiltonian is independent of the lattice parameters, and as a consequence of this, the lattice expansion coefficient:
vanishes^{3}. To assess the lattice expansion, it is thus essential to account for anharmonic effects. In this exercise, we will use the quasiharmonic approximation for this purpose ^{2}, ^{5}: In the quasiharmonic approximation, the free energy of solid is given by the total DFT energy of the electronic system and the vibrational free energy of the nuclei:
This means that our task it not to minimize the total energy with respect to the volume as you have learned it in Tutorial II by using the Murnaghan equation of state. Instead, we minimize the free energy \(F(T,V)\) with respect to volume. This needs to be done for each temperature of interest because of the temperature dependence of the vibrational free energy term \(F^\text{ha}(T, V)\). So what we need to do is:

Compute the phononic properties for slightly expanded and reduced system sizes.

find the lattice constant minimizing the free energy \(F(T, V)\) at a given temperature \(T\) using the Murnaghan equation of state.
Our procedure will be as follows: First, we generate input files for
Phonopy calculations at different volume. We will use a Python script for this called preprocess.py
, which you find in the folder:
phononswithfhivibes/Tutorial/phonons/4_QHA/input
Once the script is executed, it should create 5 working directories:
qha_35.335
qha_37.590
qha_39.939
qha_42.383
qha_44.926
In each of those you will find three input files: geometry.in
,
phonopy.in
, and aims.in
. We now need to perform a
Phonopy calculation, the Phonopy postprocess, and the reference
FHIaims calculation in each of the folders. We provide a bash script
run.sh
for this purpose that you can run with:
bash run.sh
Please, inspect the script and understand what it does. As a task: remove
every line that is not necessary to perform the actual calculations
and save as run_minimal.sh
. How many lines do you need?
Once you have computed all the thermal property files in the respective folders, we have a good dataset for vibrational free energy term \(F^\text{ha}(T, V)\).
We start to perform the postprocessing by extracting the total energies
from the aims calculations that are stored in the trajectory files
qha*/aims/trajectory.son
. The script postprocess.py
will read the trajectories in, extract the ase.Atoms
object
from there, read total energy and volume from it, and save them to a
file ev.dat
.
python postprocess.py
We can now use this information to plot a \(E\) vs \(V\) curve and fit it to
a Murnaghan equation of state. Phonopy is shipped with a set of
scripts to facilitate certain tasks. We will use the script
phonopyqha
to perform the fitting^{1}.
To fit the data, type:
phonopyqha b ev.dat eos murnaghan
This tells phonopyqha
to perform a Murnaghan fit on the data
contained in ev_total.dat
(similar to Tutorial II, where you
used the murn.py
script). The script will return the minimal
volume, energy, and bulk modulus:
phonopyqha b ev.dat eos murnaghan
> # Murnaghan EOS
> Volume: 39.959179
> Energy: 15748.028571
> Bulk modulus: 94.182465
> Parameters: 15748.028571 0.587840 4.587242 39.959179
If you run the script with the argument p
, it will also
display a plot of the energy vs. volume data points and the line
obtained from fitting the equation of state. Please try this.
We are ready to perform the final step of the quasiharmonic analysis!
Please create the directory QHA
(mkdir QHA
) and change to it (cd QHA
). From within this directory, call:
phonopyqha ../ev.dat ../qha_*/phonopy/output/thermal_properties.yaml p s
This will create a bunch of files
Cptemperature.dat
Cptemperature_polyfit.dat
Cvvolume.dat
bulk_modulustemperature.dat
dsdvtemperature.dat
entropyvolume.dat
gibbstemperature.dat
gruneisentemperature.dat
helmholtzvolume.dat
thermal_expansion.dat
volumetemperature.dat
and a plot displaying several free energy vs. volume curves at different temperatures, as well as the volume and lattice expansion coefficients vs. temperature as shown below.
Zoom into the plot displaying the \(F\) vs. \(V\) curves (helmholtzvolume.dat
) on the left panel
and take a look at the minimum of the uppermost curve marked in red. The
uppermost curve corresponds to 0 K temperature. What do you observe?
Alternatively you can compare the minimal volume you obtained from the
Murnaghan fit to the volume you find in
QHA/volumetemperature.dat
when you go to the first line (0 K).
Can you explain why the 0 K volume changes slightly when you include the
vibrational free energy?
Another interesting feature is the negative lattice expansion of Silicon below room temperature. Do you have an explanation for this? ^{4}

See https://atztogo.github.io/phonopy/qha.html#phonopyqha for reference. ↩

S. Biernacki and M. Scheffler, Phys. Rev. Lett. 63, 290 (1989). ↩

N. W. Ashcroft, N. D. Mermin, Solid State Physics, Saunders College Publishing, New York, (1976). ↩

D. S. Kim, et al., Proc. Natl. Acad. Sci. U.S.A. 115, 1992 (2018). ↩

A. Togo, L. Chaput, I. Tanaka, G. Hug, Phys. Rev. B, 81, 17430116 (2010). ↩