First slab calculation
Si(100) example calculation
Now we want to perform a simple slab calculation for the Si(100) surface. Create a new directory and create a new geometry.in
file for a 6layer slab with 40 Å of vacuum:
lattice_vector 3.8395898218429529 0.0000000000000000 0.0000000000000000
lattice_vector 0.0000000000000000 3.8395898218429529 0.0000000000000000
lattice_vector 0.0000000000000000 0.0000000000000000 46.7874999999999943
atom_frac 0.0000000000000000 0.0000000000000000 0.4274646005877638 Si
atom_frac 0.5000000000000000 0.0000000000000000 0.4564787603526582 Si
atom_frac 0.5000000000000000 0.5000000000000000 0.4854929201175527 Si
atom_frac 0.0000000000000000 0.5000000000000000 0.5145070798824473 Si
atom_frac 0.0000000000000000 0.0000000000000000 0.5435212396473417 Si
atom_frac 0.5000000000000000 0.0000000000000000 0.5725353994122361 Si
This is the same geometry that was created in Part 0: Creating slab and cluster structures for FHIaims:
Now create a control.in
file with:
# Physical settings
xc pwlda
spin none
relativistic atomic_zora scalar
# kgrid settings
k_grid 17 17 1
# Slab settings
evaluate_work_function .true.
The physical settings should look familiar to you now. We highlight some specific slabrelated settings:

k_grid 17 17 1
: The kgrid is specifically chosen for this slab system. There should be no interaction between different periodic images of the slab in z direction. Therefore, only one kpoint is needed on that axis. If you did need more than one kpoint along that axis to yield converged results, this would imply that periodic images in that direction were interacting. The kpoints along the \(k_x\) and \(k_y\) direction still have to be converged as usual. 
evaluate_work_function .true.
: This keyword turns on the calculation of the work function of the top and bottom surface. If the two sides are not equal, this should output two different values. This option requires that a reference z coordinate for the electrostatic potential evaluation deep in the vacuum is set. By default, the code will try to this set this level automatically by inspecting the geometry file. It can also be set by hand through the keywordset_vacuum_level
in thegeometry.in
file.
Now add the "light" species defaults for Si to your control.in
file via the usual command:
cat [FHIaimsdirectory]/species_defaults/defaults_2020/light/14_Si_default >> control.in
Now run the calculation (~1 minute with 4 cores) and inspect the output:
Error: The dipole correction needs vacuum zlevel to be defined
Vacuum level will be determined automatically.
Determine position of vacuum level automatically
**WARNING** : This feature is experimental and still requires some testing
**PLEASE ** : If you see it fail, please report.
Vacuum level determined:
set_vacuum_level 46.7875000000000
The code was able to set the vacuum level correctly. Always check whether this is the case for your calculation.
Once the SCF is converged, the code outputs the work function:
ESTIMATED overall HOMOLUMO gap: 0.01469701 eV between HOMO at kpoint 8 and LUMO at kpoint 72
 This appears to be an indirect band gap.
 Smallest direct gap : 0.11165629 eV for k_point 34 at 0.117647 0.882353 0.000000 (in units of recip. lattice)
However, this system has fractional occupation numbers. Since we use a finite kpoint grid,
this material is metallic within the approximate finite broadening function (occupation_type)
applied to determine the occupation numbers.
 Chemical Potential : 5.51395462 eV
Writing energy levels:
 Potential vacuum level, "upper" slab surface: 0.10357272 eV
 Potential vacuum level, "lower" slab surface: 0.10357272 eV
 Work function ("upper" slab surface) : 5.41038190 eV
 Work function ("lower" slab surface) : 5.41038190 eV
 VBM (reference: upper vacuum level) : 5.41075113 eV
 CBM (reference: upper vacuum level) : 5.39605413 eV
Selfconsistency cycle converged.
For both the top and bottom surface output we obtain the same work function with a value of \(5.41038190\) eV. This result is as expected, since we simulated a symmetric surface slab.
Work function for a semiconductor
For a semiconducting slab with a clear band gap, FHIaims arbitrarily sets the chemical potential somewhere in the band gap. The work function is calculated with respect to this chemical potential but should be calculated for a welldefined level for semiconductors, e.g. the VBM or CBM. Make sure to consider this before comparing work functions!
We can visualize this result by plotting the planaraveraged electrostatic potential. In your control.in
file, add the part:
# Potential output
output cube hartree_potential
cube origin 1.92 1.92 24.0
cube edge 39 0.1 0.0 0.0
cube edge 39 0.0 0.1 0.0
cube edge 480 0.0 0.0 0.1
cube filename electrostatic_potential.cube
This will create a cube file for the electrostatic potential called electrostatic_potential.cube
. We set the origin of the cube (roughly) in the center of the slab and then sample the cell in 0.1 Å steps, 39 steps in \(x\) and \(y\) direction and 480 in \(z\) direction. We then take the planaraverage of this potential and plot its dependence on \(z\) in the following. A postprocessing script can be found in the solution folder Tutorial/P2Slab_calculations/solutions/Si_100_workfunction
. Copy the file cube_planar_average.py
to your folder and once the cube file was created run the following python script in that folder:
cube_planar_average.py i electrostatic_potential.cube
This should create the files plavg_X.dat
, plavg_Y.dat
, and plavg_Z.dat
. You can then plot plavg_Z.dat
to visualize the planer averaged electrostatic potential in \(z\)direction.
We can also identify the results that FHIaims outputted earlier: The vaccum level is set to 0.104 eV while the chemical potential (Fermi level) is set at 5.514 eV. The work function is then defined as this difference.
The potential for the slab should in principle be symmetric. If you want to obtain a fully symmetric potential, you have to align the cube grids with the atomic centers and thus make the cube grid symmetric with respect to slab.
Remember that if we have a nonsymmetric slab, we would obtain two different work functions for the top and bottom surfaces!