Skip to content

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 file for a 6-layer 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 FHI-aims:


Now create a file with:

# Physical settings
  xc                 pw-lda
  spin               none
  relativistic       atomic_zora scalar

# k-grid 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 slab-related settings:

  • k_grid 17 17 1: The k-grid 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 k-point is needed on that axis. If you did need more than one k-point along that axis to yield converged results, this would imply that periodic images in that direction were interacting. The k-points 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 keyword set_vacuum_level in the file.

Now add the "light" species defaults for Si to your file via the usual command:

cat [FHI-aims-directory]/species_defaults/defaults_2020/light/14_Si_default >>

Now run the calculation (~1 minute with 4 cores) and inspect the output:

  Error: The dipole correction needs vacuum z-level 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 HOMO-LUMO gap:      0.01469701 eV between HOMO at k-point 8 and LUMO at k-point 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 k-point 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

      Self-consistency 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, FHI-aims 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 well-defined 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 planar-averaged electrostatic potential. In your 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 planar-average of this potential and plot its dependence on \(z\) in the following. A post-processing script can be found in the solution folder Tutorial/P2-Slab_calculations/solutions/Si_100_workfunction. Copy the file to your folder and once the cube file was created run the following python script in that folder: -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 FHI-aims 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 non-symmetric slab, we would obtain two different work functions for the top and bottom surfaces!