Surface formation energy calculation for fcc (111), (110) and (100)

Jump to: navigation, search

< Homework 1

Author(s): Laalitha Liyanage



By using the codes provided here you accept the the Mississippi State University's license agreement. Please read the agreement carefully before usage.


To calculate surface formation energy using density functional theory one should do a highly accurate bulk calculation (accurate to 1 meV) and obtain bulk energy per atom (Ebulk). Then the simulation box should be extended in the required direction to make a surface. This is equivalent to placing a vacuum over the surface. Due to periodic boundary condition two identical surfaces will be simulated by this simulation box. Therefore we only consider half of the energy required to make such a surface as the surface formation energy.

Guidelines for high accuracy bulk calculation

In VASP global precision control tag calle PREC. By setting PREC to 'Accurate' the accuracy of the calculation can be increased dramatically. Other options that control the precision in VASP include

  • ENCUT - Controls the completeness of the basis set. Could do a convergence test to determine the optimal value for your calculation. By default it will be set to the ENMAX parameter of the pseudo-potential file POTCAR.
  • KPOINT grid - Always do a k-point convergence test to determine the best k-point grid. Since the simulation box is elongated in one direction only one k-point is needed in the direction perpendicular to the surface. For other two directions k-point convergence must be done.

For more information please refer to the VASP manual[1] and VASP workshop slides.[2][3]

Guidelines for accurate surface energy calculation

The system should be relaxed (ISIF=2) with ISMEAR = 1 (Methfessel-Paxton). At the end of the relaxation run VASP will generate new positions in the CONTCAR file. Copy the CONTCAR as POSCAR. Then run a static calculation (no relaxations, NSW = 0) with the tetrahedron method (ISMEAR = -5). At the end of the static run the total energy (E0) of the system will be accurately determined. To do an accurate surface energy calculation using VASP, one has to consider several variables. To make a simulation box (POSCAR) with a surface a vacuum needs to be inserted. For convergence the following should be considered.

  • The number of atomic layer
  • The height of the vacuum
Using the following script you could generate fcc surfaces (100), (110) and (111) with your choice of the number of atomic layers and the height of the vacuum. Use the following equation to calculate surface formation energy
Surface formation eng.gif

Where Esurf is the total energy of the simulation box with the surface, N is the number of atoms in the simulation box, \epsilon is the cohesive energy per atom of the bulk structure, and A is the area of the surface.


For relaxation (geometric or ionic optimizations) always use ISMEAR = 1 and SIGMA = 0.1 or ISMEAR = 2 and SIGMA = 0.2. It is always recommended to do ionic relaxations at fixed volumes and plot the energy vs volume graph to determine the equilibrium volume. To get accurate energies obtain the structure after relaxation and run a static calculation (no relaxations NSW = 0) with ISMEAR = -5.

    • Proposed method to relax and get accurate energies
      • The system should be relaxed (ISIF=2) with ISMEAR = 1 (Methfessel-Paxton). At the end of the relaxation run VASP will generate new positions in the CONTCAR file. Copy the CONTCAR as POSCAR. Then run a static calculation (no relaxations, NSW = 0) with the tetrahedron method (ISMEAR = -5).

Example INCAR files are given below.

INCAR file for relaxation

LREAL = Auto
ENCUT = 240.3
EDIFF = 1e-6

INCAR file for static calculation

LREAL = Auto
ENCUT = 240.3
EDIFF = 1e-6

FCC surface generation script

The following python script allows you to generate super-cells with surfaces (111), (110) and (100). The input arguments are equilibrium lattice parameter of fcc bulk structure (a = lattice constant in angstroms), type of surface (surf = 100 or 110 or 111), length of the vacuum (vacuum = length in angstroms) to be inserted above the surface of the super-cell, periodicity of the super-cell (nx,ny,nz= integers), and whether or not to include an adatom (adatom = 1/0 ; true or false).

#!/usr/bin/env python
# Purpose: Calculate FCC (111), (110), and (100) surface energies
# Author: Sungho Kim and Laalitha Liyanage

import os
import sys
import math

        Usage: ./ a surf vacuum nx ny nz adatom

        Mandatory arguments
        a - equilibrium lattice constant
        surf - Type of surface 100, 110 or 111

        Optional arguments
        vacuum - length of vacuum; DEFAULT = 8.0 angstroms
        nx,ny,nz - periodicity of supercell; DEFAULT (1,1,1)
        adatom - 1/0 (True/False); DEFAULT = 0 (False)

#Default setting
vacuum = 8.0
adatom = 0

#--------------------------Surface (100)-----------------------------

def gen_data_for_fcc(a,nx=2,ny=2,nz=4):
  """ Generate datafile of FCC structure with lattice constant a """
  xa=[]; ya=[]; za=[]
  x0 = 0.0
  bx,by,bz = a*nx,a*ny,a*nz+vacuum
  x,y,z = bx,by,bz

  for i in range(nx):
    for j in range(ny):
      for k in range(nz):
        xa.append(0 + i*a); ya.append(  0 + j*a); za.append(  0 + k*a)
        xa.append(  0 + i*a); ya.append(a/2 + j*a); za.append(a/2 + k*a)
        xa.append(a/2 + i*a); ya.append(  0 + j*a); za.append(a/2 + k*a)
        xa.append(a/2 + i*a); ya.append(a/2 + j*a); za.append(  0 + k*a)
  if adatom != 0:
#        xa.append(x0); ya.append(x0); za.append(x0+nz*a)
        xa.append(bx/2.); ya.append(by/2.); za.append(x0+nz*a)
  return xa,ya,za,bx,by,bz

#--------------------------Surface (110)-----------------------------

def gen_data_for_110_fcc(a,nx=4,ny=2,nz=1):
  """ Generate datafile of FCC surface: 110:x, 112:y, 111:z """
  xa=[]; ya=[]; za=[]
  ax = a*math.sqrt(2)/2
  ay = a*math.sqrt(6)/2
  az = a*math.sqrt(3)
  x0 = 0.0
  x2 = math.sqrt(2)/4. * a
  y2 = math.sqrt(6)/4. * a
  y3 = math.sqrt(6)/6. * a
  y4 = math.sqrt(6)*5./12. * a
  y5 = math.sqrt(6)*2./6. * a
  y6 = math.sqrt(6)/12 * a
  z3 = math.sqrt(3)/3. * a
  z5 = math.sqrt(3)*2./3. * a
  bx,by,bz = ax*nx + vacuum, ay*ny, az*nz
  for i in range(nx):
    for j in range(ny):
      for k in range(nz):
        xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(x0+k*az)
        xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(x0+k*az)
        xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(z3+k*az)
        xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(z3+k*az)
        xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(z5+k*az)
        xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(z5+k*az)
  if adatom != 0:
        xa.append(x0+nx*ax); ya.append(by/2.); za.append(bz/2.)
  return xa,ya,za,bx,by,bz

#--------------------------Surface (111)-----------------------------

def gen_data_for_111_fcc(a,nx=2,ny=2,nz=4):
  """ Generate datafile of FCC surface: 110:x, 112:y, 111:z """
  xa=[]; ya=[]; za=[]
  ax = a*math.sqrt(2)/2
  ay = a*math.sqrt(6)/2
  az = a*math.sqrt(3)
  x0 = 0.0
  x2 = math.sqrt(2)/4 * a
  y2 = math.sqrt(6)/4 * a
  y3 = math.sqrt(6)/6 * a
  y4 = math.sqrt(6)*5/12 * a
  y5 = math.sqrt(6)*2/6 * a
  y6 = math.sqrt(6)/12 * a
  bx,by,bz = ax*nx, ay*ny, az*nz+vacuum
  for i in range(nx):
    for j in range(ny):
      layer = 0
      for k in range(nz):
        xa.append(x0+i*ax); ya.append(x0+j*ay); za.append(layer/3.0*az)
        xa.append(x2+i*ax); ya.append(y2+j*ay); za.append(layer/3.0*az); layer += 1
        xa.append(x0+i*ax); ya.append(y3+j*ay); za.append(layer/3.0*az)
        xa.append(x2+i*ax); ya.append(y4+j*ay); za.append(layer/3.0*az); layer += 1
        xa.append(x0+i*ax); ya.append(y5+j*ay); za.append(layer/3.0*az)
        xa.append(x2+i*ax); ya.append(y6+j*ay); za.append(layer/3.0*az); layer += 1
  if adatom != 0:
        xa.append(bx/2.); ya.append(by/2.); za.append(x0+nz*az)
  return xa,ya,za,bx,by,bz
#----------------------------POSCAR generation------------------------------------------------
def gen_poscar(xa,ya,za,bx,by,bz):
  fout = open("POSCAR","w")
  fout.write(" %22.16f  %22.16f  %22.16f\n"%(bx,0,0))
  fout.write(" %22.16f  %22.16f  %22.16f\n"%(0,by,0))
  fout.write(" %22.16f  %22.16f  %22.16f\n"%(0,0,bz))
#  fout.write("Selective Dynamics\n")
  for i in range(len(xa)):
    fout.write("%22.16f %22.16f %22.16f\n"%(xa[i],ya[i],za[i]))
#    fout.write("%22.16f %22.16f %22.16f F F T\n"%(xa[i],ya[i],za[i]))
  return len(xa)

#-------------------------------Main program---------------------------------------------------

if len(sys.argv) > 2:
        if len(sys.argv) == 3:
                a_latt = float(sys.argv[1])
                surf = sys.argv[2]
                if surf == '100' :
                        xa,ya,za,bx,by,bz = gen_data_for_fcc(a_latt)

                elif surf == '110':
                        xa,ya,za,bx,by,bz = gen_data_for_110_fcc(a_latt)

                elif surf == '111':
                        xa,ya,za,bx,by,bz = gen_data_for_111_fcc(a_latt)

        elif len(sys.argv) == 8:
                a_latt = float(sys.argv[1])
                surf = sys.argv[2]
                vacuum = float(sys.argv[3])
                nx = int(sys.argv[4])
                ny = int(sys.argv[5])
                nz = int(sys.argv[6])
                adatom = int(sys.argv[7])

                if surf == '100' :
                        xa,ya,za,bx,by,bz = gen_data_for_fcc(a_latt,nx,ny,nz)

                elif surf == '110':
                        xa,ya,za,bx,by,bz = gen_data_for_110_fcc(a_latt,nx,ny,nz)

                elif surf == '111':
                        xa,ya,za,bx,by,bz = gen_data_for_111_fcc(a_latt,nx,ny,nz)

    print "Error: wrong number of arguments!!!"
    print usage

Running the calculation

Once you generate a POSCAR file from the above script, make an INCAR and a KPOINTS file and copy the POTCAR file relevant to your group to the same directory.

POSCAR files should have more than 5 layers of atoms parallel to the interested surface. This will mean the number of atoms will be around 100 and to run this calculation you will need to execute the following command before the mpirun command.

Then do
ulimit -s unlimited
and execute VASP by
mpirun -np <no. of processors> <path of executable> 

To submit to cluster use the following method.

Submitting job to cluster

Two files are needed to submit to the cluster. One is the pbs command script and the other is a shell script to invoke ulimit command on all allocated processors. They are presented below.Both files should be in your work directory with the rest of the VASP input files.

Pbs command script

#PBS -N <name of output files>
#PBS -l nodes=4:ppn=4 
#PBS -l walltime=48:00:00
#PBS -q q64p48h@raptor
#PBS -mea
#PBS -r n

mpiexec -np 16 ./ script

ulimit -s unlimited
<Absolute path of VASP exe>
Make the file and make it executable by
chmod u+x

Extracting Data

Once the calculation is complete you will see some output files. CONTCAR gives the final geometry of the system. You could visualize CONTCAR using OVITO to make sure the structure and surface are stable. OSZICAR gives data of each electronic and ionic step. The final energy of the system (Etot) is in the OSZICAR file as well. To get the final energy execute the following command

tail -n1 OSZICAR|awk '{print $5}'


  • The surface energy should be converged in terms of k-point grid . Therefore try different k-point grids. Remember that the ratio of k-points should be inversely proportional to the lengths of the lattice vectors (in this case the edges of the simulation box). Since the supercell is elongated only one k-point is needed in the elongated direction.


Auto	        #header file
Monkhorst	#Style of Kpoints
 9  9  1	#Numbers
 0  0  0


Personal tools

Material Models