Workflow
Warning
The results here are not converged! Make sure to always check convergence of surface energies with respect to basis sets, kpoint 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:
We need to

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

Calculate the total energy of the bulk system and extract the energy \(e_{bulk}\) .

Calculate the total energy of the slab and extract the energy \(E^{slab}\) .

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 Å 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.
Task
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 pwlda
spin none
relativistic atomic_zora scalar
k_grid 12 12 12
relax_geometry bfgs 5e3
We will perform a relaxation of the atomic coordinates for a well converged kgrid of 12x12x12.
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 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 / HartreeFock s.c.f. calculation : 62993.167948580 eV
 Final zerobroadening corrected energy (caution  metals only) : 62993.167948580 eV
 For reference only, the value of 1 Hartree used in FHIaims 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 kgrid.
Task
Modify the kgrid for the slab appropriately. Keep in mind that we want to maintain the same kpoint density as in the bulk calculation!
Result
The appropriate kgrid 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\) Å while it is \(3.84\) Å in the slab! To get the same kpoint density as in the bulk we have to scale the kgrid by the lattice constant ratio: \(12\cdot 1.4 \approx 17\).
Modify the control.in
file and input the correct kgrid:
xc pwlda
spin none
relativistic atomic_zora scalar
k_grid 17 17 1
relax_geometry bfgs 5e3
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 / HartreeFock s.c.f. calculation : 62988.845773205 eV
 Final zerobroadening corrected energy (caution  metals only) : 62988.848141982 eV
 For reference only, the value of 1 Hartree used in FHIaims 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.
Task
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/Å\(^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 adirection. Copy the control.in
file for the previous slab calculation and modify the kgrid (we have a 2x1 supercell!):
xc pwlda
spin none
relativistic atomic_zora scalar
k_grid 9 17 1
relax_geometry bfgs 5e3
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.
Task
Calculate the surface energy for the Si(100) asymmetric dimer reconstruction and compare the results to the ideal surface you have calculated previously. Is the reconstruction more stable?
Result
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/Å\(^2\). The reconstruction is indeed much more stable than the ideal (100) surface!
References