Example system: Ir(100)-(2×1)-O

In this tutorial we go through a ViPErLEED structure analysis for an Ir(100) surface with adsorbed oxygen, creating a (\(2\times1\)) superstructure. The structure and analysis is based on Ref. 13. The experimental data and input model were provided by Lutz Hammer.

We will begin by running a reference calculation and a rough structure optimization. Following this, we move to a finer search grid and run a full-dynamic optimization of the inner potential \(V_{0\mathrm{i}}\). After achieving a satisfactory \(R\) factor, we use the ViPErLEED error calculation to estimate how much small changes of atomic positions affect the \(R\) factor. This will give margins for the statistical uncertainties of the determined parameters.

You can download the input files to follow along here.

Tip

If you are new to ViPErLEED, check out the more basic tutorial first, which explains the setup in more detail.

Introduction

We start from experimental LEED-I(V) curves which are saved in file EXPBEAMS.csv. As reference structure, we use a POSCAR file containing an Ir(100) surface that was hand made from a bulk structure. An additional adsorbed oxygen atom is placed at a bridge site, such that each surface iridium atom binds to one oxygen atom. The height above the surface was chosen to achieve a reasonable Ir–O bond length. See Fig. 6 for a visualization of the structure (using VESTA [12]). In this visualization, we also drew a plane orthogonal to the unit-cell vector \(\mathbf{c}\) at fractional position \(c=0.35\). We will use the BULK_LIKE_BELOW parameter to declare that everything below this cutoff is bulklike.

Fig. 6 POSCAR rendered in VESTA, viewed along \(\mathbf{a}\) (left), \(\mathbf{b}\) (center) and \(\mathbf{c}\) (right). Ir atoms are shown in blue, oxygen in red.

Hint

Of course there are other possible structural models, including adsorption at hollow, top, or Ir-sharing bridge sites. Even surface reconstructions (e.g., missing row) are conceivable. Here, we concentrate exclusively on the correct structural model, which is known from the literature [13, 14].

Often the adsorption site (or generally, the structural model) may not be known a priori. In that case, at least a rough analysis has to be performed for each reasonable model in order to find the most promising one(s), before further optimization is performed.

PARAMETERS and rough DISPLACEMENTS

We start by setting up a rudimentary PARAMETERS file, very similar to the one described in the tutorial for the Ag(100)-(1×1) surface, as per Listing 6.

Listing 6 PARAMETERS file for Ir(100)-(2×1)-O.

! ####### GLOBAL PARAMETERS #######

RUN = 1-3
THEO_ENERGIES =   50 700 3

V0_IMAG = 5.0
LMAX = 8-14

! ####### PARAMETERS FOR INTERPRETING THE POSCAR #######

BULK_LIKE_BELOW = 0.35

SITE_DEF Ir = surf top(2)
SITE_DEF O = ads top(1)

! ####### PARAMETERS FOR VIBROCC #######

T_DEBYE = 420
T_EXPERIMENT = 100
VIBR_AMP_SCALE = *surf 1.3

GLOBAL PARAMETERS

At the top of the file, we specify some general settings, such as the energy range to be used and our initial guess for the imaginary part of the inner potential \(V_{0\mathrm{i}}\). We will run a reference calculation, delta-amplitude, and a structure search back to back, so we specify RUN = 1-3.

To reduce computation time for this example, we also limit the maximum angular momentum quantum number to be used in the calculations by TensErLEED with the LMAX parameter.

PARAMETERS FOR INTERPRETING THE POSCAR

As mentioned in the Introduction, we let ViPErLEED find the bulk-repeat unit by itself. For this, we set the BULK_LIKE_BELOW parameter at the height where we drew the plane in Fig. 6.

Using the SITE_DEF parameter, we specify that the topmost iridium atoms should be treated separately. These Ir atoms are differently coordinated than bulk Ir atoms and thus will presumably have a different vibrational amplitude. For completeness, we also specify the site type O_ads for the adsorbed oxygen atom. However, since there is only one oxygen atom in the structure, this will not change the behavior and we could also skip that line.

PARAMETERS FOR VIBROCC

Since we don’t have a VIBROCC file yet, we need to specify the parameters T_DEBYE, T_EXPERIMENT, and VIBR_AMP_SCALE. Note that we can only give one value for the Debye temperature for both oxygen and iridium. The initial vibrational amplitudes will be calculated as explained in the section on T_DEBYE.

Tip

Grouping the settings in the PARAMETERS file as described above is not required, but it helps with readability.

For the delta-amplitudes calculation and structure optimization, we also need to set up a DISPLACEMENTS file. Since this is the first run starting from a totally unrelaxed model, we begin with a rather large range and a relatively rough grid, as in Listing 7.

Listing 7 DISPLACEMENTS file for the first optimization of Ir(100)-(2×1)-O.

! Label for the 1st search run
== Search z
    = GEO_DELTA
        O_ads     z = -0.20 0.20 0.01
        Ir_surf   z = -0.20 0.20 0.01
        Ir L(2-5) z = -0.20 0.20 0.01

! Label for the 2nd search run
== Search x
    = GEO_DELTA
        Ir L(1-4) xy[1 0] = -0.20 0.20 0.01

! Label for the 3rd search run
== Search u
    = VIB_DELTA
        O_ads   = -0.10 0.10 0.01
        Ir_surf = -0.03 0.03 0.01

By setting multiple search blocks in DISPLACEMENTS, we can tell ViPErLEED to execute them one after the other. When starting to optimize a new system such as here, it is generally a good idea to begin with a geometric optimization perpendicular to the surface (i.e., along \(z\)). This is because the LEED-I(V) curves are most sensitive to out-of-plane displacements.

We then follow up with rough optimizations of in-plane positions and vibrational amplitudes of surface atoms.[1] For the in-plane optimization, we use a simplified assignment for all four layers, although any movement along the \(x\) direction (here parallel to the \(\mathbf{a}\) unit vector) is forbidden by symmetry for Ir atoms in the second and fourth layers.[2] ViPErLEED automatically sorts out any symmetry-forbidden displacements, cf. Fig. 8.

Note

TensErLEED cannot simultaneously optimize displacements in different directions for a given atom. Vibrational amplitudes can be optimized together with one geometric displacement; here we do it separately to speed up the calculation. However, it is generally advisable to combine vibrational-amplitude displacements and displacements along the \(z\) direction, as the parameters may be coupled.

With the files set up, we can start the ViPErLEED calculation. During the first initialization ViPErLEED will generate IVBEAMS and PHASESHIFTS files.

Note

You will notice that the first time we run a new system, ViPErLEED stops execution after the initialization. This is on purpose, and is supposed to give the user a chance to double check the recognized symmetry and annotated POSCAR.

You will need to restart the run manually after making these checks.

When the run is finished, we will see in the log file that the \(R\) factor has decreased quite a bit. The first reference calculation gave a value of \(R_\mathrm{P} \approx 0.55\), but \(R_\mathrm{P}\) has dropped to around \(\approx 0.21\) over the search. We can get a better idea of how the search has converged by taking a look at the file Search-progress.pdf in the SUPP directory.

Fig. 7 Upper half of page 1 for Search-progress.pdf for the rough initial structure optimization of Ir(100)-(2×1)-O.

Fig. 7 shows a plot of the decreasing \(R\) factor during the various stages of the search. We clearly see that both geometric optimizations made quick progress in terms of convergence. This is expected when starting from an unrelaxed surface slab like the one used here. The other pages of the file give us some insight into how each atomic parameter developed during the structure optimization.

Fig. 8 Pages 2–4 of Search-progress.pdf for the rough initial structure optimization of Ir(100)-(2×1)-O (white space removed).

Fig. 8 shows that the (\(z\)) position of the oxygen adsorbate has changed a good amount. Further, we see that the vibrational amplitude of the oxygen atom has gone down, while the amplitude for iridium has gone up.

Important

We want to accept the optimized positions as the new starting configuration, so we need to replace our old POSCAR and VIBROCC files with the optimized ones. This can be done either manually or automatically by calling the bookkeeper utility with the --cont flag:

$ python3 bookkeeper.py --cont #[or ./bookkeeper --cont]

Fine DISPLACEMENTS and tensor-LEED error

The rough optimization has already significantly brought down the \(R\) factor. We should now continue with a finer search grid. For this, we use a similar DISPLACEMENTS file, but with much smaller range and step size (see Listing 8). We now rerun with RUN = 1-3 to perform a fresh reference calculation for the new starting positions. This is advisable because parameter deviations during the previous fit were not negligible.

Listing 8 DISPLACEMENTS file for the fine optimization of Ir(100)-(2×1)-O.

== Search z
    = GEO_DELTA
        O_def     z = -0.04 0.04 0.004
        Ir_surf   z = -0.04 0.04 0.004
        Ir L(2-5) z = -0.04 0.04 0.004

== Search x
    = GEO_DELTA
        Ir L(1-4) xy[1 0]  = -0.04 0.04 0.004

== Search u
    = VIB_DELTA
        O_def   = -0.02 0.02 0.004
        Ir_surf = -0.02 0.02 0.004

When looking at the log file after the reference calculation, we further notice something important: The \(R\) factor of the reference calculation (\(R_\mathrm{P} \approx 0.18\)) is different from the one we obtained from the superposition calculation at the end of the previous run (\(R_\mathrm{P} \approx 0.21\)). This comes from the error of the tensor-LEED approximation used for the structure optimization. In this case, the real \(R\) factor — as obtained from the reference calculation — is lower. However, this is not always the case. You should never rely on the \(R\) factor produced by the superposition calculation as a final result, but rather run a final reference calculation at the end of your analysis.[3]

Full-dynamic optimization

After the finer search run finishes, we see that the \(R\) factor has again dropped quite significantly. The \(R\) factor is now below 0.1, which already indicates very good agreement, but we can get better yet. However, before proceeding, we should accept the new best-fit structure by calling the bookkeeper utility with the --cont flag:

$ python3 bookkeeper.py --cont #[or ./bookkeeper --cont]

Now, remember that in the PARAMETERS file of Listing 6 we had to put in an initial guess for the imaginary part of the inner potential \(V_{0\mathrm{i}}\). We would now like to also optimize this non-structural parameter. However, it is not accessible in the tensor-LEED approach, which can only treat perturbations on an atom-by-atom basis. Instead, we can use the full-dynamic optimization to find an optimal value for \(V_{0\mathrm{i}}\).

During the full-dynamic optimization, multiple reference calculations are run while the chosen global parameter is varied. ViPErLEED then tries to determine an optimal value using a parabola fit. We select \(V_{0\mathrm{i}}\) for optimization by adding the following line to PARAMETERS:

OPTIMIZE V0i = step 0.5

We then chose to run the full-dynamic optimization by setting the parameter RUN = 6 and restarting. Once finished, the log file will let us know of the optimized value for the chosen parameter. ViPErLEED also produces a file called FD_Optimization.pdf in the OUT directory which contains a visualization of the \(R\) factors calculated for the various values of the parameter and a fit parabola. It is shown in Fig. 9.

Fig. 9 FD_Optimization.pdf file resulting from the full-dynamic optimization of the imaginary part of the inner potential, \(V_{0\mathrm{i}}\), for Ir(100)-(2×1)-O.

Note

ViPErLEED will also automatically add the new, optimized value to the PARAMETERS file and comment out the line containing the previous value.

Following the \(V_{0\mathrm{i}}\) optimization we can also run a final structure optimization, since the new value for \(V_{0\mathrm{i}}\) may have slightly affected the optimal positions. Using a (very fine) \(0.002\) Å grid (DISPLACEMENTS_very_fine in the provided example files) we manage to get a final \(R\) factor of around \(R_\mathrm{P} \approx 0.088\).

Error calculation

Now that we have found a good structure fit, we can run a ViPErLEED error calculation to estimate how sensitive the \(R\) factor is to small changes of specific parameters. As input for the error calculation, we need a DISPLACEMENTS file containing the desired range and steps. The format of the file is the same as used for the delta-amplitudes calculation and structure search.

Listing 9 DISPLACEMENTS file for the error calculation.

= GEO_DELTA
O_def     z = -0.05 0.05 0.002
Ir_surf   z = -0.05 0.05 0.002
Ir L(3)   z = -0.05 0.05 0.002
Ir L(4-6) z = -0.05 0.05 0.002


= GEO_DELTA                     
Ir L(2-6) xy[1 0]  = -0.10 0.10 0.004

Using the DISPLACEMENTS file in Listing 9, we run the error calculation by selecting the segment RUN = 5. The result will again be saved in the OUT directory. ViPErLEED generates a plot of the error curves in Errors.pdf and stores the raw data in Errors.csv. The Errors.pdf file in Fig. 10 shows that displacement of surface atoms even by a few picometres drastically increases the \(R\) factor. Here, atom 1 is the oxygen adsorbate and atoms 2 and 3 are the topmost iridium ones.

Fig. 10 Part of Errors.pdf showing the effects of geometric displacements.

Note

The point of intersection between the error curve for a parameter and the horizontal line labeled \(R_\mathrm{P} + \textrm{var}(R_\mathrm{P})\) gives a measure for the statistical error of the parameter [5] (see also the section on Errors).

Changes in the vibrational amplitude of the surface atoms (Fig. 11) also strongly affect the \(R\) factor.

Fig. 11 Part of Errors.pdf showing the effects of changes in the vibrational amplitude of the topmost atoms.

In general, error plots for geometric displacement tend to show a parabolic profile close to the minimum. Error plots for vibrational amplitudes tend to be more asymmetric, as these amplitudes enter the calculation differently (i.e., via the Debye–Waller factor).

Warning

Error curves are also subject to errors of the tensor-LEED approximation. \(R\)-factor values obtained for large deviations should be taken with care.

[1]

It is not recommended to fit any vibrational amplitude of bulk atoms as long as reliable information about the substrate Debye temperature is available.

[2]

Actually, no in-plane movement is allowed for these Ir atoms. They coincide with twofold rotation axes of the pmm group. See also the FreeDir column in the POSCAR file after initialization.

[3]

There is one exception to this rule. The full-dynamic reference calculation cannot provide exact results when an atom has mixed chemical composition and the elements have different optimized positions. This is because only one position can be specified for each atom. In this case, the tensor-LEED approximation is the only viable alternative.