Example system: Ag(100)-(1×1)

In this section, we go through setting up and executing a simple ViPErLEED calculation for the Ag(100)-(1×1) surface.[1]

This tutorial is designed such that you can follow along on your own system.

You can download the input files to follow along here

Tip

If this is your first time running ViPErLEED, first make sure to follow the installation instructions. For details on how to execute ViPErLEED on your system, see the How to run section.

Introduction

Our starting point for this example is:

• We have experimental LEED-I(V) curves for our surface taken at normal incidence, saved in file EXPBEAMS.csv.

• We would like to compare these intensities with a Ag(100)-(1×1) surface that was handmade from the bulk crystal structure. It is stored as a POSCAR file (see Listing 2).

In this example, we want to set up a ViPErLEED calculation and structure optimization based on these inputs alone. Our goal is to get the lowest possible \(R\) factor starting from this reference structure. If we manage to get a good agreement between the experimental and calculated curves, we can be confident that our structure model is correct.

We will start off by setting up our input files and running an initialization to check for any errors. Then, we run a reference calculation to get some initial calculated results. Afterwards, we perform a delta-amplitudes calculation and structure search to refine the results.

Initialization run and generation of phase shifts

To start off, we place the files EXPBEAMS.csv and POSCAR in our source directory. Both files have defined formats.

Listing 2 POSCAR file containing the initial structure.

Ag(100) 1x1                                     
  2.8800
    1.0000000000000000    0.0000000000000000    0.0000000000000000
    0.0000000000000000    1.0000000000000000    0.0000000000000000
    0.0000000000000000    0.0000000000000000    7.0710678100000000
   Ag  
    6  
Direct
  0.0000  0.0000  0.8000 
  0.5000  0.5000  0.7000 
  0.0000  0.0000  0.6000 
  0.5000  0.5000  0.5000 
  0.0000  0.0000  0.4000 
  0.5000  0.5000  0.3000 

We then create a PARAMETERS file (cf. Listing 3) which tells ViPErLEED what to do.

Listing 3 PARAMETERS file for the initial reference calculation.

! Comments can start with either ! or #

! ####### GLOBAL PARAMETERS #######
RUN = 0
THEO_ENERGIES = 45 700 3

! ####### PARAMETERS FOR INTERPRETING POSCAR #######
BULK_LIKE_BELOW = 0.45

SITE_DEF Ag = surf 1

! ####### PARAMETERS FOR VIBROCC #######
T_DEBYE = 330
T_EXPERIMENT = 100
VIBR_AMP_SCALE = *surf 1.3

GLOBAL PARAMETERS

Since we only want to run the Initialization for now, we set RUN = 0. For this example, we also limit our energy range to 45–700 eV, with an energy step of 3 eV (THEO_ENERGIES).

PARAMETERS FOR INTERPRETING POSCAR

We then tell ViPErLEED how to interpret the structure given in the POSCAR file. Using the BULK_LIKE_BELOW parameter, we specify that below 0.45 (unit-cell fraction along \(\mathbf{c}\)) the given structure is bulklike.

With the SITE_DEF command, we further define that the first atom in the POSCAR file (here the topmost atom) should be treated as a distinct species. See the page on the SITE_DEF parameter for details on how this works, and see also the notes on element names.

Hint

Instead of SITE_DEF Ag = surf 1, we could also set SITE_DEF Ag = surf top(1) to select the topmost silver atom, irrespective of the order in the POSCAR.

PARAMETERS FOR VIBROCC

In addition to the atomic positions, the calculation of scattering intensities also requires vibrational amplitudes for every atom in the unit cell. While the atomic positions are contained in the POSCAR file, the vibrational amplitudes are given in the VIBROCC file. However, rather than writing the VIBROCC file ourselves, we can also let ViPErLEED calculate bulk vibrational amplitudes by providing the T_DEBYE and T_EXPERIMENT parameters. By further setting VIBR_AMP_SCALE, we guess the vibrational amplitudes for the surface atoms. The VIBR_AMP_SCALE assignment means that all atoms defined as surf (via SITE_DEF) have vibration amplitudes 1.3 times larger than those calculated from the Debye temperature. All other atoms, which are implicitly assigned a def tag, will be given vibrational amplitudes determined from the Debye temperature.

That’s all the input we need to start the initialization run.

Once the run finishes, we can have a look at the log file to see if everything went as expected. Unless there was some configuration error, the log should now contain lines like these

...
Found unit cell type: square
Starting symmetry search...
Found plane group: p4m
Checking bulk unit cell...
Found SUPERLATTICE = (1x1)
...

As expected, ViPErLEED recognized our surface to be of p4m symmetry with a simple \((1 \times 1)\) termination.

During the initialization, ViPErLEED also automatically calculated electron-scattering phase shifts (based on atomic species and positions) to be used as input for the following calculations. They are stored in the PHASESHIFTS file that was copied into the source directory. This format, however, is hard to interpret for a human reader. Instead we can look at a plotted version of the same data in the file Phaseshifts_plots.pdf in the SUPP subfolder. The first page (see Fig. 3) shows the energy-dependent phase shifts for the surface atom.

Fig. 3 First page of the Phaseshifts_plots.pdf file for Ag(100)-(1×1).

The second page shows the (rather similar) phase shifts for lower-lying bulk atoms.

Reference calculation and R factor

In this simple case, we don’t need any further settings to run the reference calculation. We can just invoke viperleed.calc again after setting RUN = 1 to select the reference calculation section.

Note that the initialization is automatically executed at the start of every ViPErLEED run. Similarly, if an EXPBEAMS.csv files is provided as is the case here, the \(R\)-factor calculation is automatically executed after each reference calculation. By default, ViPErLEED uses the Pendry \(R\) factor.

Once the reference calculation finishes (only takes about 1 min with the chosen settings) we find a result for the \(R\) factor at the very end of the log file:

...
Total elapsed time: 50.78 seconds

Executed segments: 0 1 11
Final R (refcalc): 0.1732

Additionally, in the OUT directory, we find a file THEOBEAMS.csv, which contains the calculated I(V) curves and a file Rfactor_plots_refcalc.pdf, in which the experimental and calculated beams are plotted. An extract is shown in Fig. 4.

Fig. 4 First page of the Rfactor_plots_refcalc.pdf file, showing experimental (orange) and calculated (blue) I(V) curves for the first reference-calculation run for Ag(100)-(1×1).

The I(V) curves clearly show a good qualitative agreement, but the \(R\) factor of \(R_\mathrm{P} \approx 0.17\) is not great for such a simple system. We therefore proceed to the delta-amplitudes calculation and the structure search.

Note

The reference calculation also produces the tensor files which are saved in the Tensors directory. They are required as starting point for the delta-amplitude calculation and will be recognized automatically by ViPErLEED.

Delta amplitudes and structure search

To improve our \(R\) factor, we can run a local structure optimization using the tensor-LEED approach. To do this in ViPErLEED, we run a delta-amplitude calculation followed by a structure search.

First, however, we need to provide instructions about which parameters to vary in the optimization. In ViPErLEED, we give this information in the DISPLACEMENTS file (see Listing 4).

Listing 4 DISPLACEMENTS over a rough grid for Ag(100)-(1×1).

! Label for the search run
== SEARCH z

= GEO_DELTA
Ag L(1-4)   z = -0.2 0.2 0.01

See the page on the DISPLACEMENTS file for details on the syntax. Here, we allow the \(z\) positions of all silver atoms in the first four layers to vary by up to \(\pm 0.20\) Å with a step width of \(0.01\) Å.

After setting up the DISPLACEMENTS file, we can run the delta-amplitudes calculation and structure optimization by setting RUN = 2-3 in PARAMETERS. For a large system, this step can take many hours to finish, but for our simple system it only takes about 4 min to converge (using 10 CPU cores). At the end, when we take another look at our log file, we find that the \(R\) factor dropped significantly from \(R_\mathrm{P} \approx 0.17\) to \(R_\mathrm{P} \approx 0.095\). That’s not bad, but we can do a bit better.

Now that we found a better configuration, we can use the bookkeeper utility with the --cont flag to keep the new configuration and use it as our new starting point, overwriting the old POSCAR and VIBROCC files:

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

Starting from this configuration, let’s optimize with a finer grid. We change the DISPLACEMENTS accordingly, as in Listing 5

Listing 5 DISPLACEMENTS with a finer grid

== SEARCH z+u
= GEO_DELTA
Ag L(1-4) z = -0.02 0.02 0.004

= VIB_DELTA
Ag_surf = -0.02 0.02 0.004

Here we allow a \(\pm 0.020\) Å variation on a \(0.004\) Å grid. Additionally, we also allow the topmost atom to change its vibrational amplitude. This may not seem like much, but already gives \(11^5\) grid points (11 values for 5 varied parameters) and will take about three times as long as the the last run.

Important

Because we changed our reference structure, it is advisable to rerun starting with the reference calculation by setting RUN = 1-3 1. Note that we also add a second reference calculation at the end. This will remove errors due to the tensor-LEED approximation from the final result.

Once finished, we get an \(R\) factor of \(R_\mathrm{P} = 0.0836\). To visualize how our optimization went, we can also take a look at the Search-progress.pdf file in the OUT directory. The first page is shown in Fig. 5.

Fig. 5 First page of the Search-progress.pdf file showing convergence of the structure search (top) and development of atomic parameters (bottom).

The final result of the optimization (POSCAR with coordinates and VIBROCC with vibrational amplitudes) are saved in the OUT directory. Interatomic distances and angles can easily be measured in a visualization software such as VESTA [12].

Next steps

Further optimizing the structure is possible, but not very instructive. Instead we conclude this example by mentioning two other options on how to proceed with the analysis.

For a more complicated system, it may not be clear which structure parameters are most important. In this case, we could run an error calculation that shows how much the displacement of an individual atom impacts the \(R\) factor. The error calculation also provides us with an estimate for the statistical accuracy of the determined parameter values, for details see Errors.pdf.

Alternatively, we could turn to a full-dynamic optimization to also tackle parameters that are not accessible under the tensor-LEED approximation such as \(V_{0\mathrm{i}}\), the unit-cell dimensions (not relevant here, as they are well known for Ag), or the incidence angle of the electron beam.

[1]

The unpublished data was provided by Lutz Hammer.