Perfect crystal

We begin the tutorial with a simple exercise: to compute the equilibrium volume of aluminium at zero pressure and zero temperature. This requires finding the minimum of $E$ w.r.t $V$, and so our task here will be to compute the energy $E$ as function of volume.

The aluminium crystal has the face centred cubic (FCC) crystal structure. It is a cubic lattice with points at the corners of the cube and at the centres of the six faces of the cube. A possible choice of the primitive lattice vectors in this case would be ${\bf a}_1 = a(0.5, 0.0, 0.5)$, ${\bf a}_2 = a(0.5, 0.5, 0.0)$ and ${\bf a}_3 = a(0.0, 0.5, 0.5)$, i.e. three vectors pointing to the centres of three of the six faces of the cube.

We now need a computer code to calculate the energy of the crystal as function of volume. The DFT code that we will use is vasp ^9.1. vasp requires four input files, named POSCAR, INCAR, KPOINTS and POTCAR. The file POSCAR is used to define the geometry of the system, and has the following form:

a comment
-16.0  # If positive lattice parameter (Å), if negative volume of the simulation cell (Å$^3$)
0.5 0.0 0.5   # First lattice vector
0.5 0.5 0.0   # Second lattice vector
0.0 0.5 0.5   # Third lattice vector
1               # Number of atoms in the cell
Direct        # Atomic coordinates expressed in units of lattice vectors
0.0 0.0 0.0   # Atomic coordinates

The file INCAR contains several input parameters that control the accuracy of the calculations. vasp has defaults for all these parameters (for example the plane wave cutoff ENCUT, which determines the number of basis functions used in the calculations), and these default values are good enough for present purposes, so we are not going to modify any of them here beyond the couple of parameters that are already provided in the INCAR file distributed with these notes. The file KPOINTS is used to specify the grid of ${\bf k}$-points. A grid with a large number of points (high density) will give very accurate results, but it will also make the calculations more expensive. We need to find a good compromise between accuracy and cost, and increase the number of k-points until we are happy with the accuracy that we are obtaining. The KPOINTS file contains the following:

auto        # Do not change
0            # Do not change
G   # Centres the grid at $\Gamma$, makes it easier to make equivalent grids for different supercells
4 4 4       # Makes a 4x4x4 grid of points.
0 0 0       # Do not change

The file POTCAR contains the (pseudo)potential, and in this particular case the potential generated by the aluminium (pseudo)atom. Using a pseudo-potential instead of the full potential is computationally advantageous, as the inner electrons of the atoms (the core electrons) are kept frozen, because they are tightly bound to their nuclei and do not participate to chemical bonding. For this reason, it is usually a good approximation not to explicitly consider them in the calculations.

We are now ready to start exploring the energy landscape by doing a few calculations with different values of volumes: set a volume in the POSCAR file and then type:

$>$ vasp

The code will run and the energy will be available once it finishes. As we mentioned above, to obtain the equilibrium volume you can compute the energy as function of volume and look for the volume at which the energy is at its minimum. This process is facilitated by fitting the data to an appropriate function $E(V)$. This function could be a polynomial, or even a simple parabola if the calculated data only span volumes close to the minimum. A convenient function is a Birch-Murnaghan equation of state ^9.2:

$$E(V) = E_0 + \quad \quad \quad \quad \quad \quad \quad \quad \qua... ...d \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad$$

$$\frac{3}{2}V_0 K_0 \left [ \frac{3+ 6\chi}{4}\left(\frac{V_0}{V}\... ... 3\chi}{2}\left(\frac{V_0}{V}\right )^\frac{2}{3} + \frac{2\chi + 3}{4}\right],$$ (9.1)

where $E_0$ is the value of the energy at the minimum volume $V_0$, $K_0 = -V_0\left(\frac{dP}{dV}\right )_{V_0}$ the zero pressure bulk modulus and $\chi = \frac{3}{4}(4 - K_0')$, with $K_0' = \left (\frac{d K}{dP}\right )_{P=0}$. This function is defined in the file 'birch.g', and to fit it to the data we can use gnuplot. Type:

$>$ gnuplot

gnuplot$>$ lo 'birch.g'

gnuplot$>$ fit birch(x) 'filename' via emin,v0,k0

where filename is the file containing the data and emin $= E_0$, v0 $= V_0$ and k0 $= K_0$. For simplicity we do not include variations of $\chi$ in the fitting procedure, keeping $\chi = 0$.

Now repeat the study with a different number of k-points and see how your equilibrium volume changes.