Slab calculations in FHI-aims

Things to consider before running a slab calculation in FHI-aims

In this section, we want to explain how slab systems can be simulated in FHI-aims. We have already seen how slab geometries can be prepared in the previous section. Here, we summarize the most important keywords and considerations when calculating surface properties using FHI-aims. Surface will be simulated as slabs in a supercell. The slab model can be visualized as:

Before running a calculation, the key things to consider are:

• The surface is parallel to the x-y plane.

This is important for the calculation of the work function and also the dipole correction. More on this point later.

Your lattice vectors should ideally have a form like:

lattice_vector a1 a2 0.0
lattice_vector b1 b2 0.0
lattice_vector 0. 0. c3

Note that this is not always possible if the bulk system has low symmetry.

When you want to run a unit cell relaxation for a slab, it is necessary to constrain specific parts of the lattice vectors to retain the correct slab symmetry:

lattice_vector a1 a2 0.0
    constrain_relaxation z
lattice_vector b1 b2 0.0
    constrain_relaxation z 
lattice_vector 0.0 0.0 c3
    constrain_relaxation .true.

Here we constrain the z-coordinate of the first lattice vectors and all coordinates of the third one. It is crucial to constrain the last lattice vector. The unit cell relaxation would, otherwise, collapse the vacuum and severely slow down the relaxation.

When you start a FHI-aims calculation for a slab, the code will automatically detect the orientation of the slab. As an alternative to the setting lattice constraints, you can also set relax_unit_cell slab in your control.in file. The code will then automatically turn on the constraints, which is also printed in the output file:

Slab Identifier: This is Z-oriented slab and maximum thickness of the vacuum is: 40.000 Angstrom.
   Relaxation constraint has been set for slab relaxation:
   -- lattice vectors 1 and 2: Z coordinate fixed.
   -- lattice vector 3: All components fixed.

• The slab has sufficient vacuum.

Sufficient vacuum in this context means that neighboring surfaces are not interacting with each other. In many DFT codes (especially codes implementing plane-wave basis sets), you would need to run several calculations with different vacuum thicknesses to find the smallest value which gives physical results. This is because, for many DFT codes, increasing the size of the computational cell substantially increases the basis set size and thus the runtime of a calculation, even if no additional atoms are added to the system. However, in codes implementing localized basis functions (such as FHI-aims), there is negligible computational cost in adding additional vacuum to the system, so you may easily use a relatively large vacuum thickness of 50 Å or more from the onset without a noticeable performance impact.

• Charge and spin initialization.

In case of convergence problems of the SCF cycle, it is worth to think about the initialization of the partial charges of the surface atoms. Chemical intuition may help at this point (e.g. ask questions like: How many bonding partners are missing at the surface and how many electrons does the bonding partner take/donate?). The initial charges of the atoms can be assigned in the geometry.in file, like this:

atom 0.0 0.0 0.0 Mg
    initial_charge 2

The same holds for the initial_moment of an atom. For more details, we refer to the tutorial FHI-aims for Transition Metal Oxides, where this all is described in detail.

• Decide whether both surfaces, top and bottom, or only a single surface are of interest. Decouple the two surfaces electrostatically if necessary.

For most cases the system of physical interest is essentially an isolated surface, that is, there are enough “bulk” layers between the two surfaces of the slab that they are isolated from one another and are independent experimentally. If you are interested in the properties of the top surface and would like to minimize the impact of the bottom surface in your calculation, you can passivate the bottom surface and prevent physical states from appearing in the band gap. This is typically done via hydrogen or pseudo-hydrogen passivation. If the top and bottom surfaces are different, this can create a potential surface dipole in the direction normal to the surface. To compensate this dipole, you can use a dipole correction to remove this dipole ¹. In FHI-aims you can activate the dipole correction with the keyword use_dipole_correction in your control.in file.

• When relaxing the slab, should some atoms be fixed?

In certain cases it might be necessary to fix specific atoms in their position before relaxing slab, e.g. simulating the adsorption of a molecule on a surface. This can be done directly in the geometry.in file via the keyword constrain_relaxation. An example looks like this:

atom_frac 0.0000000000000000 0.0000000000000000 0.4274646005877638 Si
    constrain_relaxation .true.
atom_frac 0.5000000000000000 0.0000000000000000 0.4564787603526582 Si
    constrain_relaxation .true.
atom_frac 0.5000000000000000 0.5000000000000000 0.4854929201175527 Si
    constrain_relaxation .true.

---

1. Neugebauer, J., & Scheffler, M. "Adsorbate-substrate and adsorbate-adsorbate interactions of Na and K adlayers on Al(111)", Physical Review B, 46(24) (1992) ↩