exercise 07 :: Calculation of equilibrium structures and deformation of crystals

Next step in complexity of electronic structure calculations is the structure of the crystalline solid. Crystal samples contain replicas of the basic crystal unit cell. As typical samples used in experiments are of 1 mm size and a unit cell has typical size of 10 Angstrom the samples contain billions of replicas. To model such a number of replicas it is advantageous to consider an infinite repetition of the unit cell. This allows to solve Kohn-Sham equations within the one unit cell instead of dealing with the billions of replicas by imposing periodic boundary conditions. The idea is to expand Kohn-Sham wavefunctions and the electron density using Fourier series.

Let us consider relatively simple example, a silicon crystal. At ambient conditions it crystalizes into the diamond structure, two interpenetrating face-centered cubic (fcc) primitive lattices. The primitive unit cell contains two atoms. Further crystal information check here.

As the primitive lattice vectors and the atomic coordinates depend on the lattice parameter a, the total potential energy will be a function of the lattice parameter only (if we do not consider changes of the crystal structure) . In the crystalline solids it is convenient to introduce the cohesive energy of the solid as follows

where M is the total number of atoms in the crystal, is the total potential energy per atom (Si), and the last term is the total energy of an isolated Si atom. The cohesive energy quantifies how stable the silicon crystal is relative to the collection of isolated (non-interacting) Si atoms. The cohesive energy should go to zero for large a.

TASK 01

Calculate cohesive energy for the silicon crystal in diamond structure. For lattice parameter a consider variation from 3 to 12 Angstrom. For calculations use k-points sampling 4 4 4 and consider. Perform the calculations for three different exchange-correlation functional - local density approximation (LDA), generalized gradient approximation (GGA) by Perdew-Burke-Ernzerhof (PBE), and revised Perdew-Burke-Ernzerhof GGA (PBESOL). For pseudopotentials use the projector augmented wave (PAW) pseudopotentials:

LDA (Si.pz-n-kjpaw_psl.0.1.UPF)
PBE (Si.pbe-n-kjpaw_psl.0.1.UPF)
PBESOL (Si.pbesol-n-kjpaw_psl.0.1.UPF)

An automated bash script for self-consistent field calculations reads as follows:

#!/bin/bash          

function gener_input()
{
cat << THEEND
&control
   calculation = 'scf'
   restart_mode = 'from_scratch',
   prefix = 'silicon',
   pseudo_dir = './',
   outdir = './tmp'
 /
&system
   ibrav = 2, celldm(1) = $1, nat = 2, ntyp = 1,
   ecutwfc = 38.0, ecutrho = 152,
 /
&electrons
   diagonalization = 'david'
   conv_thr = 1.0d-8
 /
ATOMIC_SPECIES
   Si  28.086  $2
ATOMIC_POSITIONS alat
   Si  0.00  0.00  0.00
   Si  0.25  0.25  0.25
K_POINTS (automatic)
   4 4 4  0 0 0
THEEND
}

for extype in LDA PBE PBESOL
    do
      case $extype in
	LDA)		    
	      pseudo='Si.pz-n-kjpaw_psl.0.1.UPF'
	      ;;
      	PBE)
	      pseudo='Si.pbe-n-kjpaw_psl.0.1.UPF'
	      ;;
      	PBESOL)
	      pseudo='Si.pbesol-n-kjpaw_psl.0.1.UPF'
	      ;;
      esac

    echo "$extype calculations:"

    for lattconst in 3.0 4.0 4.5 5.0 5.5 6.0 6.5 7.0 8.0 9.0 10.0 11.0 12.0
       do
	  celldm=`echo "scale=6; ($lattconst)/0.5291772" | bc -l` 
          gener_input $celldm $pseudo > pw-scf.in
	  pw.x -input pw-scf.in > out.pw-scf_"$extype"_"$lattconst"
	  echo "DONE $lattconst"
       done
    done

Find lattice constants minima and cohesive energies at the minima for the three exchange-correlation functionals. Write the results to a table as shown below. Plot the cohesive energy dependencies and discuss the results.

Si (diamond)	LDA	PBE	PBESOL
cohesive energy at minimum	? eV	? eV	? eV
lattice constant at minimum	? A	? A	? A

TASK 02

Calculate cohesive energy for the silicon crystal in (beta-tin) phase using the PBESOL functional. The beta-tin phase is stable phase of silicon at high pressure (above 10 GPa). The beta-tin phase can be obtained shrinking the vertical axis by a factor c/a=0.5516. Further structural data can be found here or here. Plot and compare cohesive energies for beta-tin and diamond phases, and discuss the crystal stability.

TASK 03

To study structural equilibrium of the crystals we use enhtalpy. It allows us to predict pressure-induced phase transitions using data obtained at zero pressure and low temperatures. The enthalpy is defined as h = U + pV, where p is the pressure and V is the volume. The pressure term pV describes mechanical couplig of the crystal to the environment. At a given pressure p the most stable structure is the one with lowest enthalpy. Using results from TASK 02 calculate enthalpy as a function of volume per atom for both the diamond and beta-tin phases for pressures 0.01 GPa, 8 GPa, and 11 GPa. For conversion units of pressure use . Find pressure when both the phases are equally stable.

TASK 04

Consider a deformation of the silicon crystal studied in TASK 01 and calculate elastic constants , , and considering PBE for exchange-correlation functional. Hint: the calculations require isotropic, tetragonal and trigonal deformations. By using of Birch–Murnaghan isothermal equation of state determine bulk modulus and derivative of the bulk modulus with respect to pressure.

---

Martin Gmitra :: martin.gmitra@upjs.sk
This work is licensed under a Creative Commons Attribution 4.0 International License.