# Workflow

Warning

The results here are not converged! Make sure to always check convergence of surface energies with respect to basis sets, k-point density, number of slab layers and supercell size (1x1, 2x2, etc.).

We want to illustrate the workflow to calculate the surface energy of a stoichometric surface in detail. As shown in the introduction, we can express the surface energy for a stoichometric and symmetric slab as:

$\gamma = \frac{1}{2A} \left( E^{slab} - N\cdot e_{bulk} \right) \, .$

We need to

1. Create the bulk structure. Also create the symmetric (if possible) slabs for the desired surface and extract the surface area $$A$$.

2. Calculate the total energy of the bulk system and extract the energy $$e_{bulk}$$ .

3. Calculate the total energy of the slab and extract the energy $$E^{slab}$$ .

4. Calculate the surface energy $$\gamma$$.

This workflow looks rather simple but there are some potential pitfalls that we will highlight when they can come up.

We want to illustrate each point for the Si(100) surface.

## Si(100) example

### Step 1: Creating the slab and bulk structures

Create a directory called Si_100_surface_energy and enter it. In there we create 2 folders: bulk and slab. You should now have the following folder structure:

Si_100_surface_energy
├── bulk
├── slab


Enter the bulk folder and create the Si bulk structure geometry.in file with:

lattice_vector 5.4299999999999997 0.0000000000000000 0.0000000000000000
lattice_vector 0.0000000000000000 5.4299999999999997 0.0000000000000000
lattice_vector 0.0000000000000000 0.0000000000000000 5.4299999999999997
atom_frac 0.0000000000000000 0.0000000000000000 0.5000000000000000 Si
atom_frac 0.7500000000000000 0.7500000000000000 0.7500000000000000 Si
atom_frac 0.0000000000000000 0.5000000000000000 0.0000000000000000 Si
atom_frac 0.7500000000000000 0.2500000000000000 0.2500000000000000 Si
atom_frac 0.5000000000000000 0.0000000000000000 0.0000000000000000 Si
atom_frac 0.2500000000000000 0.7500000000000000 0.2500000000000000 Si
atom_frac 0.5000000000000000 0.5000000000000000 0.5000000000000000 Si
atom_frac 0.2500000000000000 0.2500000000000000 0.7500000000000000 Si


Now go back and enter the slab folder and create the geometry.in for the Si(100) 4 bilayer slab with 40 Å of vacuum:

lattice_vector 3.8395898218429529 0.0000000000000000 0.0000000000000002
lattice_vector -0.0000000000000002 3.8395898218429529 0.0000000000000002
lattice_vector 0.0000000000000000 0.0000000000000000 43.4399999999999977
atom_frac 0.5000000000000000 0.5000000000000000 0.5156250000000000 Si
atom_frac 0.5000000000000000 0.0000000000000000 0.5468750000000000 Si
atom_frac 0.0000000000000000 0.0000000000000000 0.5781250000000000 Si
atom_frac 0.0000000000000000 0.5000000000000000 0.6093750000000000 Si
atom_frac 0.5000000000000000 0.5000000000000000 0.3906250000000000 Si
atom_frac 0.5000000000000000 0.0000000000000000 0.4218749999999999 Si
atom_frac 0.0000000000000000 0.0000000000000000 0.4531250000000000 Si
atom_frac 0.0000000000000000 0.5000000000000000 0.4843749999999999 Si


Visualize both structures (e.g. with GIMS) and make sure that they look sensible. The structures were created with the pymatgen library and a similar script as the one that was used in P0 (This should link to P0 somehow). Now calculate the surface area for Si(100) with $$A = \left| \vec{a} \times \vec{b} \right| \, ,$$ where $$\vec{a}$$ and $$\vec{b}$$ are the lateral lattice vectors of the slab.

Calculate the surface area for the Si(100) slab.

Result

You should obtain $$A\approx 14.74$$$$^2$$.

### Step 2: Calculation for the bulk system

We now want to start the calculation for the bulk system. Go back in the bulk folder and also create the corresponding control.in file:

xc                 pw-lda
spin               none
relativistic       atomic_zora scalar

k_grid   12  12  12

relax_geometry bfgs 5e-3


We will perform a relaxation of the atomic coordinates for a well converged k-grid of 12x12x12. Now add the "light" species defaults for Si to your control.in file via the usual command:

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


Now run the calculation (~1 minutes with 4 cores). It should only take a single relaxation step. Once the calculation is finished, your output should resemble:

------------------------------------------------------------------------------
Final output of selected total energy values:

The following output summarizes some interesting total energy values
at the end of a run (AFTER all relaxation, molecular dynamics, etc.).

| Total energy of the DFT / Hartree-Fock s.c.f. calculation      :         -62993.167948580 eV
| Final zero-broadening corrected energy (caution - metals only) :         -62993.167948580 eV
| For reference only, the value of 1 Hartree used in FHI-aims is :             27.211384500 eV
| For reference only, the overall average (free atom contribution
| + realspace contribution) of the electrostatic potential after
| s.c.f. convergence is                                          :            -12.847000415 eV


The total energy of the bulk is $$E_{bulk}=-62993.167948580$$ eV. What we want to have for the surface free energy is the total energy per formula unit. The Si bulk structure has 8 atoms, so our desired total energy per formula unit is $$e_{bulk}=E_{bulk}/8$$.

### Step 3: Calculation for the slab system

Now we need to do the calculation for the slab system.

Warning

Because we are calculating energy differences, it is crucial that we calculate all structures at the same level of theory with the same computational settings!

Copy the control.in file from the bulk folder to the slab folder. We still need to modify one computational setting in the slab case: the k-grid.

Modify the k-grid for the slab appropriately. Keep in mind that we want to maintain the same k-point density as in the bulk calculation!

Result

The appropriate k-grid in x and y direction for the slab needs to be scaled according to the lattice vectors. Notice that the lattice constant in the bulk is $$5.43$$ Å while it is $$3.84$$ Å in the slab! To get the same k-point density as in the bulk we have to scale the k-grid by the lattice constant ratio: $$12\cdot 1.4 \approx 17$$.

Modify the control.in file and input the correct k-grid:

xc                 pw-lda
spin               none
relativistic       atomic_zora scalar

k_grid   17 17 1

relax_geometry bfgs 5e-3


Run the calculation for the slab (~2 minutes for 4 cores). Your output should resemble:

------------------------------------------------------------------------------
Final output of selected total energy values:

The following output summarizes some interesting total energy values
at the end of a run (AFTER all relaxation, molecular dynamics, etc.).

| Total energy of the DFT / Hartree-Fock s.c.f. calculation      :         -62988.845773205 eV
| Final zero-broadening corrected energy (caution - metals only) :         -62988.848141982 eV
| For reference only, the value of 1 Hartree used in FHI-aims is :             27.211384500 eV
| For reference only, the overall average (free atom contribution
| + realspace contribution) of the electrostatic potential after
| s.c.f. convergence is                                          :             -3.210854747 eV


We can extract the total energy of the slab as $$E_{slab}=-62988.845773205$$ eV.

### Step 4: Calculate the surface energy

We now have all the ingredients to calculate the surface energy.

Calculate the surface energy for the Si(100) 4 bilayer slab using the obtained results.

Result

The formula for the surface energy in this case is: $$\gamma = \frac{1}{2A} \left( E^{slab} - 8\cdot e_{bulk} \right) \, ,$$ since we have 8 atoms in the slab. You should obtain $$\gamma\approx0.147$$ eV/Å$$^2$$.

## Si(100) reconstruction

The ideal Si(100) surface that we have just simulated is actually a metastable structure. The stable Si(100) structure is an asymmetric dimer reconstruction [1]. We now want to test whether this reconstruction actually has a lower surface energy than the ideal surface. Since the reconstruction occurs in a 2x1 supercell, we have to create a 2x1 supercell of the ideal Si(100) surface. This can be done with the structure builder in GIMS after you have inserted the ideal surface geometry as shown here:

Only relaxing the ideal surface doesn't lead to the reconstruction. We have to slightly disturb the top and bottom most Si atoms to get out of the metastable state. Create a new folder called dimer_reconstruction and create a new geometry.in file with:

lattice_vector 7.679180 0.000000 0.000000
lattice_vector -0.000000 3.839590 0.000000
lattice_vector 0.000000 0.000000 43.440000

atom_frac 0.250000 0.500000 0.515625 Si
atom_frac 0.250000 0.000000 0.546875 Si
atom_frac 0.000000 0.000000 0.578125 Si
atom_frac -0.020000 0.500000 0.609375 Si
atom_frac 0.270000 0.500000 0.390625 Si
atom_frac 0.250000 0.000000 0.421875 Si
atom_frac 0.000000 0.000000 0.453125 Si
atom_frac -0.000000 0.500000 0.484375 Si
atom_frac 0.750000 0.500000 0.515625 Si
atom_frac 0.750000 0.000000 0.546875 Si
atom_frac 0.500000 0.000000 0.578125 Si
atom_frac 0.520000 0.500000 0.609375 Si
atom_frac 0.730000 0.500000 0.390625 Si
atom_frac 0.750000 0.000000 0.421875 Si
atom_frac 0.500000 0.000000 0.453125 Si
atom_frac 0.500000 0.500000 0.484375 Si


Notice how the top and bottom most atoms are symmetrically displaced by 0.02 in a-direction. Copy the control.in file for the previous slab calculation and modify the k-grid (we have a 2x1 supercell!):

xc                 pw-lda
spin               none
relativistic       atomic_zora scalar

k_grid   9 17 1

relax_geometry bfgs 5e-3


Now you can run the calculation but note that this is more computationally demanding (~15 minutes with 40 cores). You can also download the output from the git solutions folder. Once the calculation is done, you can visualize the relaxation path with GIMS:

We can see that it starts with a symmetric dimer structure and then relaxes to the asymmetric dimer reconstruction.

The formula for the surface energy in this case is: $$\gamma = \frac{1}{4A} \left( E^{slab} - 16\cdot e_{bulk} \right) \, ,$$ since we have 16 atoms in the slab and a 2x1 supercell. You should obtain $$\gamma\approx0.09$$ eV/Å$$^2$$. The reconstruction is indeed much more stable than the ideal (100) surface!