Tuesday Lab – June 5, 2001
Electronic Structure of atoms, molecules and solids using VASP
Instructor – Blair Tuttle (PSU – Erie)
I. Atoms:
1. Hydrogen atom
ENCUT = 136 eV 
LDA 
GGA 
Notes 
TOTEN (eV) 
0.88 
1.10 
ISPIN = 1, raises E by ~0.9 eV for LDA; ENCUT à up, TOTEN à more neg by ~ 0.05 eV 
TOTEN + EATOM (eV) 
13.01 
13.53 
2. Silicon atom:
ENCUT = 136 eV 
LDA 
GGA 

TOTEN (eV) 
0.70 
0.85 

Notes 
EATOM = 115.7 
EATOM = 116.73 
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.781 
0.766 
0.74 
K(LDA) = 36.7 eV/Ang^2 K(GGA) = 35.9 eV/Ang^2 
TOTEN min (eV) 
6.585 
6.639 
===== 

w (cm^{1}) 
4460 
4420 
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.28 
2.20 
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.05 
< 0.0 unbound 
LDATOTEN (H3) = 7.52 eV 
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) 
18.79 

TOTEN – SiH_{3} (eV) 
13.43 

E_{B} (eV) Exp = 3.95 eV 
4.22 
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. TOTEN (LDAISPIN=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 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) 
5.57 
4.57 
4.63 
Band Gap (eV) 
0.49 
0.64 
1.12 
Notes: TOTEN GAMMA 0.8 X LDA 11.88 5.03 5.52 GGA 10.88 5.06 5.60 
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) 
0.812 

Binding energy (eV) 
2.1 

w (cm^{1}) 
4230 
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?
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 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) 
197.86 

TOTEN – H^{0} (eV) 
192.46 

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}) = 0.4 eV
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?
Yes, neutral hydrogen is a shallow donor ie the defect level is near the conduction band minimum. At the gamma kpoint, 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” 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