Skip to content

Creating slabs with ASE or Pymatgen

Before simulating any kind of system, we have to create an appropriate starting structure. This is usually doable for (known) bulk structures. The problem can be much more difficult in surface systems. Cutting a certain surface direction from an arbitrary bulk system is a non-trivial task. Luckily, there are several libraries/utilities available that can create surface structures with little input. In this part we want to highlight some convenient ways to create slab structures for use in FHI-aims.

Warning

You have to be careful when creating surfaces from any starting structure in an automized way. Miller indices are typically defined via the conventional lattice vectors. If you use the general surface utilities listed below, make sure to start from the conventional bulk structure to avoid any discrepancies with available literature.

Creating surface slabs using ASE

Creating surface slabs using ASE is rather intuitive. We will focus on the example of the Si(100) surface. Here is an example python script to create a Si(100) surface using ASE's built-in utility functions:

from ase.build import diamond100

# Create the diamond 100 surface of Si directly using an utility function
sur = diamond100("Si", [1,1,6], a=5.43, vacuum = 20.)

sur.write("Si_100.in", scaled=True)

The first line imports the necessary function from ASE.

from ase.build import diamond100

The diamond100 utility function directly constructs the desired Si(100) surface. The ASE website provides also additional predefined surface functions and instructions for building arbitrary surfaces..

We can then directly construct the desired Si(100) surface in a single line:

# Create the diamond 100 surface of Si directly using an utility function
sur = diamond100("Si", size=[1,1,7], a=5.43, vacuum = 20.)

We create a Si(100) slab with 1x1x7 times the size of the minimal unit cell and a vacuum of 20 Å in z-direction on both sides of the slab. The total amount of vacuum in this case is 40 Å!

The structure can then directly be written to file in FHI-aims format:

sur.write("Si_100.in", scaled=True)

Running the above script should create the file Si_100.in which can be visualized, e.g. in VESTA:

vesta_Si_100

Creating surface slabs using pymatgen

You can also use the pymatgen library to construct surface slabs. The implemented surface utilities are slightly more generalized but also more complex compared to the ones in ASE. When creating surfaces with a specified Miller index, we can generally get several terminations that describe the same surface direction but have a different surface composition. Pymatgen can output all symmetrically inequivalent terminations simultaneously.

Here, we only show a simple example. More details can be found at pymatgen surface page and the mentioned references.

We create the same Si(100) surface as before. Notice that we start from the conventional bulk structure. An example script is:

from pymatgen.core.surface import SlabGenerator, Lattice, Structure
from pymatgen.io.ase import AseAtomsAdaptor

# Create the diamond Si bulk first
lattice = Lattice.cubic(5.43)
Si = Structure(lattice, ["Si", "Si", "Si", "Si",
                         "Si", "Si", "Si", "Si"],
               [[0.00000, 0.00000, 0.50000],
                [0.75000, 0.75000, 0.75000],
                [0.00000, 0.50000, 0.00000],
                [0.75000, 0.25000, 0.25000],
                [0.50000, 0.00000, 0.00000],
                [0.25000, 0.75000, 0.25000],
                [0.50000, 0.50000, 0.50000],
                [0.25000, 0.25000, 0.75000]])

# Now we use the pymatgen SlabGenerator to create the Si(100) surface 
slabgen = SlabGenerator(Si, (1,0,0), min_slab_size=2, min_vacuum_size=6, in_unit_planes=True, center_slab=True)
Si_100 = slabgen.get_slabs()[0]

# Pymatgen doesn't support the aims input/output format yet
# We use ase to do the geometry output
ase_py = AseAtomsAdaptor()
ase_struc = ase_py.get_atoms(Si_100)

ase_struc.write("Si_100.in", scaled=True)

As you can see, there is slightly more effort involved to create the same surface slab. Note, however, that the surface utilities in pymatgen are formulated for arbitrary crystal symmetry. This also leads to the slightly higher degree of complexity.

We start by creating the Si diamond structure by hand.

# Create the diamond Si bulk first
lattice = Lattice.cubic(5.43)
Si = Structure(lattice, ["Si", "Si", "Si", "Si",
                         "Si", "Si", "Si", "Si"],
               [[0.00000, 0.00000, 0.50000],
                [0.75000, 0.75000, 0.75000],
                [0.00000, 0.50000, 0.00000],
                [0.75000, 0.25000, 0.25000],
                [0.50000, 0.00000, 0.00000],
                [0.25000, 0.75000, 0.25000],
                [0.50000, 0.50000, 0.50000],
                [0.25000, 0.25000, 0.75000]])

The SlabGenerator class then creates the surface structures based on the input structure and additional keywords.

slabgen = SlabGenerator(Si, (1,0,0), min_slab_size=2, min_vacuum_size=6, in_unit_planes=True, center_slab=True)
Si_100 = slabgen.get_slabs()[0]

The second input is the desired Miller index, (100) in our case. Refer to the pymatgen wiki page for the additional inputs. There is only 1 surface termination in the case of Si(100), but this is not generally the case as pointed out in the introduction.

In the end, we have to convert the pymatgen Structure to an ASE Atoms object and output the slab as a FHI-aims geometry.in file.

ase_py = AseAtomsAdaptor()
ase_struc = ase_py.get_atoms(Si_100)

ase_struc.write("Si_100.in", scaled=True)