Purpose: To give the student an appreciation for the accuracy and speed of
state-of-the-art 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 state-of-the-art electronic structure code.  It uses a planewave basis set (and therefore periodic boundary conditions), ultra-soft pseudo-potentials and solves the Kohn-Sham equations using iterative techniques.  For exchange-correlation, both the LDA and GGAs are available.

The basic input files are:

POSCAR (lattice vectors and basis atoms)

KPOINTS (both user choices and Monkhorst-Pack)

POTCAR (contains psuedopotential (PP) info include LDA/GGA and atom info)

Note each sub-directory 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 sub-directories, 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 Perdew-Wang 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

1. Enter ”lab_tuttle/_h” and look at the files esp. INCAR.
2. Run executable script “>RUN” which runs vasp and writes the total energy (TOTEN) to a file called “ENERGY”
3. To compare the TOTEN result to experiment (E = -13.60 eV) you have to add the reference energy which is EATOM=  -12.13 for LDA calculations as you can find in the OUTCAR file. To perform GGA calculations set GGA flag in INCAR file.  Note EATOM = -12.43 eV for GGA (91) functional.
4. Vary ENCUT and ISPIN and see how TOTEN varies.

2.      Silicon atom:

1. enter “lab_tuttle/_si1” ; repeat above

This result will be used later when we calculate the cohesive energy for bulk c-Si.  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₂

1. Enter ”lab_tuttle/_h2” and look at the files esp. INCAR.
2. Run executable script “>RUN” which calculates the energy of H₂ for several bond lengths and write the results to the file “EvsD”
3.  View “EvsD” and record TOTEN and D(H-H) for minimum by inspection or by fitting data to a parabola.  A simple way to get a fit to data is to use “xmgr” which is available on ncsa machines.  See appendix A for instructions on fitting data with “xmgr.” From a parabolic fit you can extract a harmonic frequency as demonstrated along with unit conversion in Appendix B.
4. Test the convergence of the results with respect to ENCUT and ISPIN.

5. Calculate the binding energy per hydrogen for H₂ (see lecture notes):

1.      –E_B = ( Etot(H₂) – 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₂  the zero point energy, Ezp ~ 0.26 eV.

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.    

1. Enter ”lab_tuttle/_h3” and look at the files esp. INCAR.
2. Run executable script “>RUN” which calculates a gradient search for the minimum energy configuration within LDA, see file = ENERGY.  Note this calculation may take a while since the potential energy surface is flat. You may want to kill the job before it finishes.
3. Activate GGA flag in the INCAR to see that no new minimum will be found. In this case, there is a barrier of about 0.8 eV for the exchange of H. Although VASP does not have a transition state search implemented, in this case, the potential energy surface is so simple that you could find the barrier “by hand.”
4. Calculate the binding energy per H atom assuming Ezp ~ 0.0 eV. 

3.      Si-H molecules: SiH₄ and SiH₃ – Here we will calculate the Si-H binding energy in silane (SiH₄).

1. Enter “lab_tuttle/_sih”
2. Run executable script “>RUN” which calculates the minimum energy configurations for SiH₄ and SiH₃ in LDA. The results are in ENERGYsih4 and ENERGYsih3, respectively.     
3. Report total energies below and calculate Si-H binding energy using the Si-H zero point energy is Ezp ~ 0.24 eV.

d.   For SiH₃, see the effects of spin polarization by performing a spin-averaged 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 c-Si in its (lowest energy) FCC lattice.

1. Enter “lab_tuttle/_csi2”
2. Run executable script “>RUN” which calculates the total energy and then the bandstructure for Si (see vasp web page for details on band calculation). The results are in files ENERGY and BANDS.  In the BANDS file, kpt#1 is the L point, kpt#6 is the Gamma point and kpt#11 is the X point.
3. You can use xmgr to view BANDS and compare to the band structure found by the empirical pseudo-potentials in R. Martin’s lab.
4. Calculate the band gap and cohesive energy:

– E_coh = [E_tot(c-Si) – 2E_tot(Si) +Ezp]/2 :for E_tot(Si), see I.2.a, and Ezp ~ .05 eV.

2.      H₂ in c-Si.  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₂ in semiconductors generally.

1. Enter “lab_tuttle/_csi8h2”
2. Run executable script RUN which calculates the energy of the 8 atom c-Si cell and then calculates the energy for varying H₂ bond lengths. The output data is in files ENERGYcsi8, ENERGYcsi8H2 and  EvsD, respectively.
3. Using the EvsD data file, we can determine the minimum bond length and the minimum energy for H₂ in silicon. 
4. Calculate bond length, binding energy and frequency for H₂ in Si. To find the binding energy use:

–E_B = ( Etot(c-Si+H₂) – E(c-Si) – 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₂ in a vacuum (II.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 sub-directory “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 sub-directories for “lab_tuttle/_si32donor” In particular, look at the eigenvalues in OUTCAR or EIGENVAL.  The highest occupied eigenvalue in the H⁰ 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” sub-directories.  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 c-Si 32 cell. I have entered this value in the table below, but you can find it in the OUTCAR file in the sub-directory: lab_tuttle/_si32/_csi/_gamma

c.       From the notes the donor level is the Fermi energy (E_F^+/0) when

E_form(H⁰) = E_form(H⁺¹)

and with some algebra we can solve for

E_F^+/0 =  E_tot(H⁰) - 1 E_VBM – E_tot(H⁺¹)

d.      Compare E_F^+/0 to the LDA band gap (see III.1.c) and to the defect eigenvalue for H⁰ (see III.3.a). Is neutral hydrogen a shallow donor? If yes, then can we use neutral atomic hydrogen to dope silicon p-type?

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” sub-window 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 à Non-Linear Curve Fitting}

                        {“Non-Linear Curve Fitting” window should appear}

                        Click in “Formula” sub-window

                                    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₁ m₂ / (m₁ + m₂ )

w (cm^-1) = 521.4 SQRT( A2 [units = eV/Ang^2 ] / m [au] )  [ units = cm^-1]

1 au = 1 atomic unit 

For H₂:   m = m₁ m₂ / (m₁ + m₂ ) ~ 0.5 au