Summer School on Computational Materials Science

SOLUTIONS

Tuesday Lab – June 5, 2001

Electronic Structure of atoms, molecules and solids using VASP

Instructor – Blair Tuttle (PSU – Erie)

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.   TOTEN (LDA-ISPIN=2) = -13.68 so DE = 0.25 eV

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?

Yes, the paircorrelation function is flat after the first nearest neighbor peak.  Liquid silicon still shows some covalency which is why you have a first neighbor peak. There are only 32 atoms in the simulation which is why the function is jagged but a larger model would be quit smooth after r =3A.

d.      If you quenched the melted supercell, would it go to a crystalline or amorphous state?

I expect the model to go to a highly defective amorphous state… 

e.       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⁺¹) = 0.4 eV

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?

Yes, neutral hydrogen is a shallow donor ie the defect level is near the conduction band minimum. At the gamma k-point, the neutral defect eigenvalue is close to the donor level (there is a fair amount of scatter in the 32 atom calculations so other kpoints are not as close).  This is because there is little steric relaxation between the two configurations so most of the donor energy is associated with occupying the H⁺ defect level.

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