Purpose: To give the student an appreciation for the accuracy and speed of
stateoftheart electronic structure methods based on density functional
theory.
Tuesday Lab – June 5, 2001
Electronic Structure of atoms, molecules and solids using VASP
Instructor – Blair Tuttle (PSU – Erie)
Introduction:
In the exercises described below we will explore the properties of atoms, molecules and solids using the Vienna ab initio simulation package (VASP). You will have the opportunity to calculate bond lengths, vibrational frequencies, bond/cohesive energies, defect levels, and band structures. At any time, feel free to skip around or try your own calculations. The goal of this lab is simply to introduce you to some of the types of calculations that can be done and to give you a feeling for the convergence and the accuracy of the results.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
You need an executable version of VASP in order to do the following lab
See VASP web site: http://cms.mpi.univie.ac.at/vasp/Welcome.html
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
Getting Started:
1. First download file “lab_tuttle_cms01” from Summer School web site.
2. unzip tar file: (>gunzip lab_tuttle_cms01.tar.gz).
3. Untar directory tree: (> tar –xvf lab_tuttle_cms01.tar).
4. Now, there should be in a local directory, a directory called “lab_tuttle”
5. Move to subdirectory (>lab_tuttle/_vasp) then link your vasp executable to “vasp”
Notes on VASP:
At this point it would be useful to get a web browser to display the VASP user guide:
http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html
I will give a brief introduction to the file structure. You may want to look at the VASP guide for more details and to see what the default settings are.
VASP is a stateoftheart electronic structure code. It uses a planewave basis set (and therefore periodic boundary conditions), ultrasoft pseudopotentials and solves the KohnSham equations using iterative techniques. For exchangecorrelation, both the LDA and GGAs are available.
The basic input files are:
POSCAR (lattice vectors and basis atoms)
KPOINTS (both user choices and MonkhorstPack)
POTCAR (contains psuedopotential (PP) info include LDA/GGA and atom info)
Note each subdirectory will have a POTCAR file. These are LDA PPs.
INCAR (contains most info regarding what type of run will be performed)
The important output files are:
OUTCAR (main output file contains energy/force info for each ionic step)
CONTCAR (final coordinates and velocities; has same format as POSCAR)
EIGENVAL (eigenvalues for each kpt in reciprocol coordinates)
PCDAT (ave. pair correlation, see below for more info)
DOSCAR (density of states, you need large Kpt sampling for reasonable results)
CHGCAR (final charge density on real space grid determined by the FFT spacing; if you are fluent in Matlab, you can visualize the charge density; see my fortran program charge.f which reads the charge in and write out the charge in one direction)
In most subdirectories, I have written a script that runs the vasp executable and distills some info from the run. We will mostly concern ourselves with total energies in this lab which are reported in the OUTCAR file In my scripts I have created a file called ENERGY which contains just the total energies TOTEN from the OUTCAR file, as shown here:

alpha Z PSCENC = 0.002052
Ewald energy TEWEN = 2.042824
1/2 Hartree DENC = 6.022587
V(xc)+E(xc) XCENC = 2.386085
PAW double counting = 0.000000 0.000000
entropy T*S EENTRO = 0.000000
eigenvalues EBANDS = 7.253506
atomic energy EATOM = 12.046704

free energy TOTEN = 0.884076 eV
===============================================
The free energy TOTEN is the sum of the values above it. In the labs, we will focus on this energy and the RUN scripts you will notice often “grep TOTEN.”
See the vasp web site for a complete list of flags you can vary in the INCAR file. Here, I will mention a few flags of interest to the calculations we will do:
ICHARG = how charge is initially calculated: 1 = file from previous run, 2 = atomic orbitals (this is what we will use mostly), 11 = keep charge fixed (this is used for bandstructure and DOS calculations)
ENCUT = cutoff energy for plane wave basis set (eV) –most of the calculations I have set up are for 10 Ryd or 136 eV basis cutoff which allows the jobs to run in just a minute or so. Feel free to test the convergence of the results.
ISPIN = spin polarization? 1= No; 2 = Yes
GGA = 91 :calculates GGA (using functional of PerdewWang from 1991) energies using lda PPs (it is better to use GGA PPs but this gives a reasonable picture of how GGA differs from LDA and is less cumbersome than copying PPs files ….).
I. Atoms:
1. Hydrogen atom
ENCUT = 136 eV 
LDA 
GGA 
Notes 
TOTEN (eV) 

TOTEN + EATOM (eV) 
2. Silicon atom:
ENCUT = 136 eV 
LDA 
GGA 

TOTEN (eV) 

Notes 
This result will be used later when we calculate the cohesive energy for bulk cSi. There is no experiment to compare with here so there is no need to add EATOM to TOTEN.
II. Molecules:
1. Hydrogen diatomic molecule: H_{2}
ENCUT =136 
LDA 
GGA 
Exp. 
Notes 
D(HH) min. (Ang) 
0.74 

TOTEN min (eV) 
===== 

w (cm^{1}) 
4400 
1. –E_{B} = ( Etot(H_{2}) – 2E(H) + Ezp )/2 where
2. Ezp = Sum [(1/2) hbar w] = 6.2 x 10^{6} cm Sum[w (cm^{1})] eV
3. Note for H_{2 } the zero point energy, Ezp ~ 0.26 eV.
ENCUT =136 
LDA 
GGA 
Exp. 
Notes 
E_{B} (eV) 
2.26 
2. Hydrogen triatomic molecule: As mentioned in the lectures, LDA tends to overbind which leads to situations where overcoordination is erroneously predicted. The triatomic hydrogen molecule is one example.
ENCUT =136 
LDA 
Exp. 
Notes 
E_{B} (eV) 
< 0.0 unbound 
3. SiH molecules: SiH_{4} and SiH_{3} – Here we will calculate the SiH binding energy in silane (SiH_{4}).
ENCUT = 136 eV 
LDA 
Notes 
TOTEN – SiH_{4} (eV) 

TOTEN – SiH_{3} (eV) 

E_{B} (eV) Exp = 3.95 eV 
d. For SiH_{3}, see the effects of spin polarization by performing a spinaveraged LDA calculation. To do so, set ISPIN = 1 in the INCAR file.
III. Solids:
1. Crystal Silicon: Here we will calculate both the cohesive energy and the bandstructure for cSi in its (lowest energy) FCC lattice.
– E_{coh} = [E_{tot}(cSi) – 2E_{tot}(Si) +Ezp]/2 :for E_{tot}(Si), see I.2.a, and Ezp ~ .05 eV.
ENCUT = 136 eV 
LDA 
GGA 
Experiment 
Cohesive E (eV) 
4.63 

Band Gap (eV) 
1.12 

Notes: 
2. H_{2} in cSi. Using an 8 atom supercell, we will see the effects of a semiconductor on the properties of a molecule. The results here even with a low plane wave cut off and a small supercell are indicative of the effects seen for H_{2} in semiconductors generally.
–E_{B} = ( Etot(cSi+H_{2}) – E(cSi) – 2E(H) + Ezp )/2
To find the harmonic frequency fit data in EvsD file using xmgr (as in II.1.c).
e. Compare to your LDA calculations above for H_{2} in a vacuum (II.1).
ENCUT = 136 eV 
LDA (cSi) 
Notes 
Bond length (Ang) 

Binding energy (eV) 

w (cm^{1}) 
3. Melting Silicon: Here I have simulated the melting of silicon by annealing a 32 atom supercell at T = 4000 K for 2 ps compared to annealing silicon at T = 300 K for 1 ps.
a. Enter subdirectory “lab_tuttle/_si32melt”
b. View the pair correlation for both the melt and the crystal (>ghostview paircorr.ps) The average pair correlation is taken over the last 1 ps of the melt simulation.
c. Has the model fully melted compared to the crystal data? If you quenched the melted supercell, would it go to a crystalline or amorphous state?
d. All the files for the simulations can be found in the subdirectories.
3. Atomic hydrogen adsorbed in Si: Here we will calculate the donor level for atomic hydrogen in bulk crystalline silicon using a 32 atom super cell. Actually, you only have to look at the results. Since the calculations take between 10 – 20 minutes,I have already performed the calculations.
a. Please look at the files in each of the subdirectories for “lab_tuttle/_si32donor” In particular, look at the eigenvalues in OUTCAR or EIGENVAL. The highest occupied eigenvalue in the H^{0} cell is the defect level which is approximately the donor level.
b. To calculate donor level we need:
i) the total energy for the neutral and positive charge state. These are in the ENERGY files of “_ho” and “_hp” subdirectories. Fill in the table below.
ii) the value of the valence band maximum for the 32 atom supercells which is approximately (to within 0.05 eV) the valence band maximum (VBM) of the cSi 32 cell. I have entered this value in the table below, but you can find it in the OUTCAR file in the subdirectory: lab_tuttle/_si32/_csi/_gamma
ENCUT = 136 eV 
LDA 
Notes 
TOTEN – H^{+1} (eV) 

TOTEN – H^{0} (eV) 

E_{VBM} (eV) 
+5.00 
c. From the notes the donor level is the Fermi energy (E_{F}^{+/0}) when
E_{form}(H^{0}) = E_{form}(H^{+1})
and with some algebra we can solve for
E_{F}^{+/0}_{ }= E_{tot}(H^{0})  1 E_{VBM }– E_{tot}(H^{+1})
d. Compare E_{F}^{+/0} to the LDA band gap (see III.1.c) and to the defect eigenvalue for H^{0 }(see III.3.a). Is neutral hydrogen a shallow donor? If yes, then can we use neutral atomic hydrogen to dope silicon ptype?
Appendices:
Appendix A: Using xmgr
Here is a step by step as to how to fit a polynomial to your energy versus bond length data “file=EvsD” using xmgr:
>{enter directory with EvsD file}
>xmgr &
{ graphical xmgr window should appear}
{click on: File à Read à Sets }
{“Read Sets” window should appear}
Scroll in “File” subwindow for “EvsD”
Click on “EvsD”
Click on “Autoscale on Read”
Click on “OK”
{data should appear on back window}
Click on “Cancel”
To fit data to a parabola:
{Click on: Data à Transformations à NonLinear Curve Fitting}
{“NonLinear Curve Fitting” window should appear}
Click in “Formula” subwindow
Enter formula: y = A0 + 0.5*A2*(x – A1)^2
Set “Number of Parameters” to 3
Set A1 = 0.7 {this initial guess is needed for fitting to work well}
Click on “Run 100 Steps”
{Make sure fitting looks good; if not adjust A0, A1 and A2}
{“Results” window should appear}
A0 = minimum energy (eV)
A1 = minimum bond length (Ang)
A2 = spring constant (eV/ Ang^2)
{see below for using A2 to get w (cm^{1}) }
Appendix B: Calculating w (cm^{1}) from A2 (eV/Ang.^2)
I will skip most of the derivation which involves several unit transformations:
w (cm^{1}) = f (sec^{1}) / 2 p c
f (sec^{1}) = sqrt ( k/ m ) = sqrt ( A2/ m )
m = m_{1} m_{2} / (m_{1 }+ m_{2 })
w (cm^{1}) = 521.4 SQRT( A2 [units = eV/Ang^2 ] / m [au] ) [ units = cm^{1}]
1 au = 1 atomic unit
For H_{2}: m = m_{1} m_{2} / (m_{1 }+ m_{2 }) ~ 0.5 au