Calculating Lattice Parameters

One of the first steps in modeling a solid material is to determine the lattice parameters for the material. Solid materials have periodic building blocks, the unit cell, that make up the whole crystal structure. Correctly modeling the unit cell means we can model a material. Incorrectly modeling the unit cell means we are adding artificial strain to the system and any results we have may be suspect. Note that the unit cell in the real world, e.g. from experiment, may be different than the unit cell we calculate using density functional theory (DFT) or other methods. DFT is only an approximation of the real world, so the most stable, strain-free unit cell may be slightly different than the real-world equivalent.

The size and shape of the unit cell is described by three vectors, the lattice vectors, as well as the positions of the atoms within the unit cell. As shown below, the angles between the vectors and vector lengths are key to describing the unit cell. See [wikipedia](https://en.wikipedia.org/wiki/Crystal_structure) for more description. Usually (but not always) we have the simple case where all the angles are at 90^o. For this case we only need to find the lattice lengths, or lattice parameters.

Case I - All lattice parameters equal

This is the simplest case. An example would be CsCl, shown below. We only need to determine the length of the box to describe the unit cell. For CsCl the atoms are either in the middle of the cell or at the corners. A procedure for this is given below.

• Obtain the crystal structure, usually a cif file.

Many online crystallography databases exist with crystal structures. The cif file contains information on the lattice vectors, space group, atomic positions, etc.

• Create the input files for the crystal structure.

For VASP, this means creating the POSCAR file, as well as INCAR, KPOINTS, and POTCAR files. In CP2K you only need an input file.

• Model many different possible lattice parameters.

You may have to create different input files and run them all, like creating different POSCAR files by hand and running them all. Or you can create a script to automate this procedure. Below is a bash script for Turing that will loop over several lattice parameter values (the "for" loop). It creates a POSCAR with the lattice value, runs vasp, takes the output and appends to SUMMARY.acell, makes backups of output names, and repeats. Running this is a nice way to get different energies for different unit cell sizes.

#!/bin/bash
#SBATCH -p short 
#SBATCH -J CsCl-lattice
#SBATCH -N 1
#SBATCH --time=4:00:00
#SBATCH --mem 40000
#SBATCH --ntasks-per-node=18
#SBATCH --cpus-per-task=1
#SBATCH --gres=gpu:0
export OMP_NUM_THREADS=1 
module unload cp2k/gcc-4.8.5/2.6.0
module load vasp/5.3
cd /home/bob/work3/vasp/CsCl-lattice
rm WAVECAR
for i in 4.58 4.59 4.60 4.61 4.62 4.63 4.64 ; do
cat >POSCAR <<!
fcc:
  1.0
$i 0.0 0.0
0.0 $i 0.0
0.0 0.0 $i
  1 1
direct
0 0 0
0.5 0.5 0.5
!
echo "a= $i" 
mpirun -n 18 vasp > log-$i
E=`tail -1 OSZICAR` ; echo $i $E >>SUMMARY.acell
cp OUTCAR OUTCAR-$i
cp OSZICAR OSZICAR-$i
done
cat SUMMARY.acell

• Make a plot of energy versus the lattice parameter.

Once you have energies of different lattice parameter, you need to plot them. Excel or any plotting program works fine. You'll get a plot like the one below. The stable unit cell has the lattice parameter with lowest energy. Use this lattice parameter for future calculations. In the example below 4.13 Å is the stable lattice parameter.

• Backup your files.

Always do this. :)

Note on number after decimal! I usually like lattice parameters with accuracy to two numbers after the decimal. This is usually good enough. Sometimes getting accuracy to one number will work, like 4.1 Å. I often may start with 'rough' calculations with only one number, then refine once I know around where the minimum should be.

Make sure you have enough data. Make sure you have enough data below and above the proposed minimum so you can actually prove that the minimum is the correct minimum. For example, if 4.13 Å is the proposed minimum, we'd want at least 4.12 Å and 4.14 Å also to show that 4.13 Å is a minimum.

Case 2- Two lattice parameters equal

This case involves looking for two independent lattice parameters that minimizes the system energy. It's similar to first case, with just a little more complexity. An example is TiO₂ where there are two lattice parameters to find, a and c.

• Get the crystal structure, cif file.

• Create the input files for the crystal structure.

• Model many different possible lattice parameters.

The following script works great.

#!/bin/bash 
#SBATCH -p short
#SBATCH -J tio2-rut-lat
#SBATCH -N 1
#SBATCH --time=4:00:00
#SBATCH --mem 40000
#SBATCH --ntasks-per-node=18
#SBATCH --cpus-per-task=1
#SBATCH --gres=gpu:0
export OMP_NUM_THREADS=1
module unload cp2k/gcc-4.8.5/2.6.0
module load vasp/5.3
cd /home/bob/work/vasp/tio2-lattice
rm WAVECAR
for i in 4.58 4.59 4.60 4.61 4.62 4.63 4.64 4.65 4.66 4.67 4.68 4.69 4.70; do
for j in 2.92 2.93 2.94 2.95 2.96 2.97 2.98 2.99 3.00 3.01 3.02 3.03 3.04 3.05 ; do
cat >POSCAR <<!
tio2:
  1.0
$i 0.0 0.0
0.0 $i 0.0
0.0 0.0 $j
  2 4
direct
0 0 0
0.5 0.5 0.5
0.305 0.305 0
0.695 0.695 0
0.195 0.805 0.5
0.805 0.195 0.5
!
echo "a= $i , c = $j"
mpirun -n 18 vasp > log-$i-$j
E=`tail -1 OSZICAR` ; echo $i $j $E >>SUMMARY.acell
cp OUTCAR OUTCAR-$i-$j
cp OSZICAR OSZICAR-$i-$j
done
done
cat SUMMARY.acell
cut -d ' ' -f 1-2,5 SUMMARY.acell > summary-energies-TiO2.txt

• Make a plot of energy versus the lattice parameters.

A contour plot is a great way to show the energy versus the two lattice parameters. You can use a program like Sigmaplot or matplotlib. See the plot below. In this case the minimum energy looks to be around 4.68 and 3.03.

• Verify that you have a minimum.

Look that exact numbers to determine where the minimum energy lies. Also make sure that you have all the energies near the proposed minimum so that you can verify that you actually are at the minimum. In this case we'd want the energies for (4.67,3.03), (4.68,3.03), (4.69,3.03)(4.67,3.02)(4.68,3.02)(4.69,3.02),(4.67,3.04),(4.68,3.04),(4.69,3.04). If we have these energies we can verify by looking at the numbers that (4.68,3.03) is the actual minimum.

• Backup your files.

Case 3- More complicated cases

For cases with all three lattice parameters different or angles other than 90^o the situation is more complicated, since you're looking for minimum over 3 or more parameters. The situation, however is similar to previous cases...you simply need to get energies for different parameters and find the parameters with minimum energy. Visualization is more difficult and you need to pay closer attention to the data.