Example system: Fe₂O₃(1−102)-(1×1)

Hematite (ɑ-Fe₂O₃) is one of the most common, naturally occurring iron oxides. Its surfaces are popular model systems for (single-atom) catalysis [15, 16], photocatalysis and the study of mineral–water interfaces.

The \((1\bar{1}02)\) surface — also referred to as \((012)\), or r-cut — exhibits two stable terminations under UHV conditions: a stoichiometric, bulk-truncated, (1×1) termination, and a reduced (2×1) reconstruction [15]. Here we use the bulk-truncated (1×1) surface as an example system of an oxide surface with medium complexity and will compare experimental LEED-I(V) data with a theory calculation in ViPErLEED.

This example is a simplified version of the example calculation shown in the main ViPErLEED paper [3]. For more details on the full calculation, see the paper and the supplementary information.

All input files needed to follow along on your own machine can be found here.

Tip

If you are new to ViPErLEED, check out this tutorial first, which explains how to set up a calculation in more detail.

Note

If you want to follow along on your own machine, note that the calculation may take a long time to complete. On a modern 48-core machine, each looped structure optimization took about 6–10 hours to finish.

Introduction

To set up the calculation, we need to start from a structure model of the (1×1) surface in the POSCAR format. We can easily build this ourselves from a (ɑ-Fe₂O₃) bulk structure. For this example, we use the structure by Maslen et al. [17], which can be found on the Crystallography Open Database (entry 2101167).

Warning

Bulk structures can be found in various online databases, but you should be careful to choose a model with the correct experimental lattice parameters (as determined by XRD). Lattice parameters from, e.g., DFT calculations may be off by ~1%, which will strongly affect the \(R\) factor. Lattice parameters are not accessible in the tensor-LEED approximation, but they can be optimized in a full-dynamic optimization.

To go from the hexagonal bulk unit cell to the unit cell for the \((1\bar{1}02)\) surface, we need to apply a suitable change of basis

\[(\mathbf{a}', \mathbf{b}', \mathbf{c}') = (\mathbf{a},\mathbf{b},\mathbf{c}) P,\]

with

\[\begin{split}P = \frac{1}{3} \begin{pmatrix} 3 & 1 & 2 \\ 0 & 2 & 4 \\ 0 & -1 & 1 \end{pmatrix} .\end{split}\]

In VESTA, you can do this by selecting Edit ‣ Edit Data ‣ Unit cell ‣ Transform. To avoid problems later on, it is useful to also delete the unit cell symmetry before applying a transformation. You can do that by pressing the Remove symmetry button in Edit ‣ Edit Data ‣ Unit cell. When applying the transformation, VESTA may prompt you and ask what to do with the atoms in the new unit cell. Make sure to choose Search atoms in the new unit-cell and add them as new sites. The unit-cell transformation will introduce duplicate atoms, which can be deleted using the Remove duplicate atoms button in Edit ‣ Structure parameters.

To create a suitable surface slab (as expected by ViPErLEED), we can then replicate the unit cell along \(\mathbf{c}\) and (optionally) remove the topmost layer to create a vacuum interface. In VESTA you can do this by applying a transformation that “stretches” the cell in the \(\mathbf{c}\) direction. You may have to use the “Initialize current matrix” button first to apply the transformation to the previously transformed unit cell. Use

\[\begin{split}\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 2 \end{pmatrix}\end{split}\]

to make the slab twice as thick along \(\mathbf{c}\). Choose again Search atoms in the new unit-cell and add them as new sites. You can then add a vacuum gap by applying the same transformation again, this time not adding any new atoms (choose Do nothing in the prompt). This is not strictly required by ViPErLEED, but makes it easier to recognize which layer is our surface.

Our starting POSCAR file contains 40 atoms in total, which corresponds to four repeat units (layers) of 4 Fe and 6 O atoms each. See Fig. 12 for a side view of the structure.

PARAMETERS and VIBROCC

As in the previous examples, we begin with a simple PARAMETERS file (Listing 10).

Listing 10 PARAMETERS for ɑ-Fe₂O₃\((1\overline{1}02)\)-(1×1).

# GLOBAL PARAMETERS
RUN = 0

BEAM_INCIDENCE = THETA 0.0, PHI 90.0

# POSCAR INTERPRETATION
SITE_DEF O = surf top(2)
SITE_DEF Fe = surf top(2)

BULK_LIKE_BELOW = 0.35

# DEBYE TEMPERATURE FOR VIBRATIONAL AMPLITUDES
T_DEBYE = 480
T_EXPERIMENT = 300

The settings concerning global parameters and the VIBROCC file have been discussed in detail elsewhere [see, e.g., the Ag(100) example]. In particular, we specify the parameter BULK_LIKE_BELOW with a suitable value (0.35 for the POSCAR file constructed in the Introduction). ViPErLEED will use it to derive the values of the BULK_REPEAT, LAYER_CUTS, and N_BULK_LAYERS parameters, as shown in section Initialization below.

We also use the SITE_DEF parameter to define which sites should be treated separately by ViPErLEED. In this case, we declare two surface sites for iron and oxygen, selecting the two topmost atoms of each species. These atoms will be assigned dedicated vibrational amplitudes and electron–atom scattering phase shifts (see file PHASESHIFTS). All remaining atoms are automatically given a default label instead.

Note also that we set the BEAM_INCIDENCE parameter to \(\phi=90°\) to match the input structure with the crystal orientation used in the experiment. This setting won’t have any effect initially, since the polar angle is set to \(\theta = 0°\). However, we will optimize \(\theta\) later on, at which point the value of \(\phi\) will be important.

To set initial vibrational amplitudes for the surface atoms, we can either supply a VIBROCC file, or let ViPErLEED calculate them based on the Debye temperature of the sample. In this case, we set the Debye temperature and experiment temperature via the T_DEBYE, T_EXPERIMENT parameters. Note that this Debye temperature was chosen based on literature values and an previous optimization run. For more details, see the main ViPErLEED publication.

Initialization

We can now run the initialization to check if all input files are interpreted correctly. If so, the log should look similar to Listing 11.

Listing 11 Log file after initialization for ɑ-Fe₂O₃\((1\overline{1}02)\)-(1×1).

...
Reading structure from file POSCAR
ViPErLEED is using TensErLEED version 1.76.

STARTING SECTION: INITIALIZATION

Loaded EXPBEAMS.csv file containing 41 beams (total energy range: 1.2e+04 eV).
Found unit cell type: rectangular
Starting symmetry search...
Found plane group: pg[0 1]
Detected bulk repeat vector: [-2.51774 2.33620 3.68229]
Checking bulk unit cell...
Found SUPERLATTICE = (1x1)
Starting bulk symmetry search...
Found bulk plane group: pg[0 1]
Generating phaseshifts data... 
Generating BEAMLIST...
Finishing section at 11:55:31. Section took 11.79 seconds.

Starting cleanup...
Wrote manifest file successfully.

Finishing execution at 2024-02-20 11:55:31
Total elapsed time: 11.83 seconds

Executed segments: 0

The initialization log shows that the input files are interpreted correctly and that the right plane group (pg[0 1]) is detected. It also shows that the bulk structure and repeat vector have been identified. This is reflected in the PARAMETERS file, which should now contain additional lines as in Listing 12.

Listing 12 PARAMETERS after initialization.

...
! BULK_LIKE_BELOW = 0.35            ! line commented out automatically

! ######################################################
! #  THE FOLLOWING LINES WERE GENERATED AUTOMATICALLY  #
! ######################################################

BULK_REPEAT = [-2.51774 2.33620 3.68229] ! Automatically detected repeat vector
LAYER_CUTS = 0.1250 0.2500 0.3750
N_BULK_LAYERS = 1

The detected layers and bulk repeat vector are also shown in Fig. 12 for clarity.

Fig. 12 Visualization of the crystal structure with highlighted layers, and the bulk repeat vector. Fe atoms are depicted as large yellow spheres, O atoms as small red ones. This figure is taken from the first ViPErLEED paper [3].

At this stage, before proceeding, you should always check the detected bulk structure and spacing, as it is crucial for the calculation. To help with this, ViPErLEED provides the files POSCAR_bulk and POSCAR_bulk_appended. The former contains the crystal structure of the bulk as detected by ViPErLEED; the latter contains the surface slab with additional bulk layers appended to it. You can visualize these files in VESTA to check if the bulk structure is as you expect. In particular, one should make sure that the bulk stoichiometry for POSCAR_bulk, and the bulk layer spacing for POSCAR_bulk_appended, are correct.

Reference calculation and structure optimization

If the initialization ran without errors, we can proceed to the reference calculation and structure optimization. You can find more details about what happens in each section of a ViPErLEED LEED-I(V) calculation here.

In short, the reference calculation performs a full-dynamic LEED calculation to generate I(V) curves for the input structure. I(V) curves are very sensitive to the exact atomic positions and vibrational amplitudes. Because of this, it is necessary to optimize the structure to get the best fit to the experimental I(V) data. The reference calculation is the starting point for this optimization, but full-dynamic calculations are computationally too expensive for the optimization, which requires a large number of iterations. We can perform the optimization in a computationally more efficient way by using a perturbative approach, where we calculate the effect of small changes to the input structure. This is done in the delta-amplitudes calculation and structure optimization.

To run these sections, we also need to specify which parameters we want to vary (i.e., the parameter space). We do this by providing a DISPLACEMENTS file, as in Listing 13.

Listing 13 DISPLACEMENTS

<loop>

== SEARCH L1-2 z+vib
  = GEO_DELTA
  * L(1-2) z = -0.05 0.05 0.01

 = VIB_DELTA
  *surf = -0.02 0.02 0.005 

== SEARCH L1-2 xy
  = GEO_DELTA
  * L(1-2) xy = -0.03 0.03 0.01

</loop>

Here, we specify that we want to optimize the atomic out-of-plane positions of all atoms in the first and second layer from −0.05 Å to +0.05 Å in steps of 0.01 Å. At the same time, we also optimize the vibrational amplitudes of the surface atoms. We subsequently turn to the in-plane position of the atoms in the first two layers. Note also the <loop> and </loop> tags at the beginning and end of the file, which indicate that these steps should be repeated until convergence. For more details on the syntax, see the page on the DISPLACEMENTS file.

We can now start the calculation by setting the RUN parameter to 1-3 1 and starting viperleed.calc again. This will execute the reference calculation, the delta-amplitude calculation and the looped structure optimization in order. By appending the 1 to the 1-3, we also tell ViPErLEED to finish with a second reference calculation of the optimized structure.

Once finished, you can check the log files to see if the calculation ran without errors. You can also find the final \(R\) factor at the end of the log file (see the extract in Listing 14).

Listing 14 Final log

...
Finishing execution at 2024-02-22 00:53:04
Total elapsed time: 9:47 hours

Executed segments: 0 1 11 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 
31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 31 12 2 3 
31 12 2 3 31 12 1 11
Final R (refcalc): 0.4023
Final R (superpos): 0.4936

Fig. 13 shows the evolution of the \(R\) factor during the optimization process. The initial \(R\) factor (\(R_\mathrm{P}\approx0.78\)) is very high. This is not unusual when starting from a totally unrelaxed, bulk-truncated structure, as multiple iterations of reference calculation and structure optimizations are usually necessary for convergence.

Fig. 13 Upper half of page 1 for the Search-progress.pdf file for the first optimization run of ɑ-Fe₂O₃\((1\overline{1}02)\)-(1×1).

Note

You may also notice that the \(R\) factor for the second reference calculation is quite different from the \(R\) factor obtained at the end of the structure optimization. This is to be expected, as the perturbative tensor-LEED approximation used in the structure optimization is only valid for small changes to the structure. The final reference calculation removes this error and gives a more accurate \(R\) factor.

Optimizing the incidence angle

In our case, the \(R\) factor has improved significantly after the first structure optimization, but is still far from a good fit at \(R_\mathrm{P} \approx 0.40\). This is because our measurement was performed on a UHV manipulator that does not have enough degrees of freedom in rotation to perfectly align the sample with the LEED optics. We thus know that we need to optimize the incidence angle of the electron beam to get a closer correspondence between the experimental and calculated I(V) curves.

The incidence angle is not a parameter that is accessible in the tensor-LEED approximation. In ViPErLEED, we can instead optimize the incidence angle using a full-dynamic optimization, which optimizes a parameter by performing a full-dynamic reference calculation at each iteration.

To do this, we need to edit the PARAMETERS file (see Listing 15) to set RUN=6 and provide the OPTIMIZE parameter to specify which parameter we want to optimize.

Listing 15 PARAMETERS, with OPTIMIZE specified for the polar incidence angle \(\theta\) of the primary electron beam.

RUN = 6

BEAM_INCIDENCE = THETA 0.0, PHI 90.0

OPTIMIZE theta = step 1, minpoints 5, convergence 0.1
...

Here, we set the initial step size to 1° and convergence criterion to 0.1°. After running the calculation, we find a minimum at about \(\theta \approx 0.87°\) which already reduces the \(R\) factor to about \(R_\mathrm{P} \approx 0.22\). We can also see that the I(V) curves, and consequently the \(R\) factor, are very sensitive to the incidence angle in the FD_Optimization.pdf and FD_Optimization_beams.pdf files (see Fig. 14 and Fig. 15).

Fig. 14 \(R\) factor as a function of the polar incidence angle as shown in FD_Optimization.pdf.

Fig. 15 I(V) curves for beam (1|1) as shown in FD_Optimization_beams.pdf.

Note

If the incidence angle is off-normal, it is advisable to optimize it relatively early on, before performing structural optimization deeper into the surface. Otherwise, the error in the I(V) curves may incorrectly be compensated by moving the atoms relative to the incident beam, which will result in incorrect positions relative to the bulk. This becomes more likely if bulklike layers contribute only weakly to the I(V) curves, i.e., when optimizing layers far below the surface.

Refining the structure

After optimizing the incidence angle, we can perform another structure optimization to refine the structure. To do this, we again provide a DISPLACEMENTS file to specify the range of parameters we want to optimize. We can then run the calculation by setting the RUN parameter back to 1-3 1. Note also that we should now remove or comment out the OPTIMIZE parameter, as the superfluous parameter will raise a warning otherwise.

At this point, the best strategy to achieve the optimal fit is generally not known a priori, and the process may involve some trial and error. In particular, you may need to play around with multiple iterations of coarse- and fine-grained structure optimizations, and possibly also adjust the convergence parameters.

In this case, we found that by first performing another looped structure optimization for the first two layers (atoms * L1-2), and then adding the third-layer atoms (atoms * L1-3) to the optimization, we could reliably reduce the \(R\) factor to about \(R_\mathrm{P} \approx 0.16\).[1]

For reference, the final \(R\) factor obtained by in the more detailed calculation discussed in the main ViPErLEED paper was \(R_\mathrm{P} = 0.1537\).

[1]

Optimizing all three layers right away with the default convergence parameters caused the \(R\) factor to get trapped in a local minimum.