ViPErLEED package I: Calculation of $I(V)$ curves and structural optimization

Low-energy electron diffraction (LEED) is a versatile, quick and easy-to-use characterization technique routinely used in many surface-science laboratories. Acquiring LEED images typically takes minutes, limited more by the time it takes to move a sample into position than the actual measurement. The resulting LEED pattern offers a reciprocal-space fingerprint of the surface unit cell, where each “spot” corresponds to a beam of diffracted electrons that fulfills the constructive-interference Bragg condition of the 2D periodicity of the surface. These data can indicate whether a given sample preparation yields the desired surface structure. However, this snapshot of the reciprocal lattice is only a fraction of the full information available from LEED. Spot-profile analysis (SPA)LEED [1, 2], as well as analysis of the background [3], can offer information about defects and surface order. LEED $I(V)$ can yield accurate positions for the atoms in the unit cell. In LEED $I(V)$, the intensities $I$ of the scattered electron beams are measured as a function of electron energy, tuned via the accelerating voltage $V$. Since these “$I(V)$ curves” can also be calculated theoretically, experimental LEED-$I(V)$ data can be directly compared to a structural model of a surface. Comparison between calculations and experiment is performed through the so-called $R$ factor, which evaluates agreement across the entire beam set and breaks it down into a single value [4, 5]. Thus, with a relatively quick and simple laboratory-based experiment, one can confirm or reject hypothetical surface models, perform local optimizations of the atomic coordinates in such a model, or obtain a robust fingerprint to confirm the reproducibility of sample preparation [6, 7].

Crucially, the hardware for these experiments is already present in many surface-science laboratories, and only minor adjustments are required to upgrade most standard LEED setups to LEED-$I(V)$ capability. The bigger hurdles are the analysis of the data and the calculation of the theoretical $I(V)$ curves based on a given model. ViPErLEED, presented here, is a comprehensive package that combines a simple hardware add-on with software for acquiring and extracting the $I(V)$ data, and for calculating and comparing to theoretical $I(V)$ curves. Figure 1 shows the workflow with input and output.

The whole package consists of three parts, which are presented in three separate publications: The data acquisition (hardware), discussed in Part III [8]; the data extraction (spot tracking) and processing, discussed in Part II [9]; and the LEED-$I(V)$ calculation and optimization, presented here. Specifically, this work discusses the software to calculate $I(V)$ curves for a given reference structure and to optimize this structure based on the match with experimental data. While these three parts are optimized to be used together, they all employ standard file formats (such as comma-separated values — CSV) for input and output and can be used independently.

An extensive treatment of the theory of quantitative LEED, as well as the standard computational treatment, can be found in Refs. 10, 4, 11, 12, 13, 14, 7, 15, 16, 17, 6. These are only outlined briefly here to introduce important terms and concepts used in this work.

LEED is highly surface sensitive due to the low inelastic mean free path of electrons in solids ($\approx$5–15 Å in the relevant energy range of up to 1000 eV [18]). Because of the strong electron–solid interaction, electrons typically experience multiple scattering events before leaving the surface [19]. Thus, calculations of LEED-$I(V)$ spectra need to include paths with multiple scattering events. Crucially, small displacements (on the order of a few picometres) strongly affect the intensities of scattered beams. The high sensitivity to atomic positions gives LEED $I(V)$ the potential for accurate determination of atomic coordinates. This, however, can only be achieved through a careful optimization process while comparing calculated and experimental $I(V)$ curves.

Scattering of low-energy ($\approx 30$–2000 eV) electrons with one atom in a solid is well approximated as occurring in a spherically symmetric potential (muffin-tin approximation) and conserving the angular momentum $\ell$. Single scattering events of electrons with energy $E$ at a given site can then be described by partial-wave phase shifts $\delta_{\ell}(E)$. Electron–phonon scattering is included as an effective Debye–Waller factor, which depends on the vibration amplitudes of the atoms. Inelastic scattering is treated via an imaginary part of the inner potential, leading to an exponential decay with distance. The real part of the inner potential is defined by an approximate function, derived together with the scattering phase shifts with a boundary condition to the vacuum [25, 26]. With these few approximations, one can perform full-dynamic calculations of the multiple scattering of electrons, yielding theoretical $I(V)$ curves that can be compared to experiment [10].

In a typical LEED calculation, the surface is split into “layers”, which should be separated along the surface normal by distances of at least $1$ Å. Each layer can consist of one or multiple atoms. Multiple scattering within each layer is then evaluated in angular-momentum space, resulting in layer-diffraction matrices. Propagation and scattering between the layers can be treated in a plane-wave representation [14, 10]. Figure 2 shows a typical layer definition for an example system of medium complexity: the $\alpha\text{-Fe}_{2}\text{O}_{3}(1\overline{1}02)$-$(1\hskip 1.00006pt\times\hskip 1.00006pt1)$ surface [24]. This system will be referred to several times below, and it is one of the examples given in the ViPErLEED documentation (see the Supplemental Material [27]). This layering scheme also allows for straightforward treatment of scattering far below the surface: One or two layers are defined as a bulk repeat unit (e.g., Layer 4 in Fig. 2) with a given bulk repeat vector. This bulk unit can then be stacked through repeated doubling until the norm of the reflection and transmission matrices of the bulk remains unchanged between successive steps [14].

The “full-dynamic” LEED calculations discussed so far are computationally demanding. The most expensive part is the inversion of matrices needed for intralayer scattering, which scales as $d^{n}$ ($n=2.38$–3, depending on the algorithm) for $d\times{}d$ matrices. Intralayer scattering matrices have $d=N(\ell_{\text{max}}+1)^{2}$ for $N$ atoms in the layer and a maximum angular momentum $\ell_{\text{max}}$ [7]. As discussed in Section I.2, the high sensitivity to atomic coordinates on the picometre scale implies that basic comparisons of theoretical models to experimental data require local optimization of these coordinates, entailing numerous calculations. The number of optimization parameters also scales with system size, with a concomitant increase in the number of trial calculations [28].

The so-called tensor-LEED approximation makes these optimizations more efficient: A full-dynamic calculation is performed only for a reference structure; changes to the diffraction intensities through small modifications of the reference structure are treated by first-order perturbation theory [11, 12, 13]. This approach is ideal for structural optimization aimed at improving the agreement between theoretical and experimental $I(V)$ curves, as the full-dynamic calculation only needs to be performed once. The optimization can then be performed using much cheaper evaluations of a set of perturbations, which scale linearly with the number of parameters under variation [12]. Importantly, the initial structural guess must be close to the target structure: The tensor-LEED approximation breaks down for displacements of more than $\approx 0.4$ Å [12], and significant errors are already introduced much sooner. This means that new reference calculations must be performed when changes to the original structure become too large. It is worth noting that the tensor-LEED approximation is not limited to purely geometrical perturbations of the reference structure. Changes of atomic vibration amplitudes [29] and statistical chemical substitutions [30] (including vacancies) can be treated with the same method and even concurrently [7, p. 135].

Several implementations exist for calculating LEED-$I(V)$ curves, most of which are derived from the pioneering works of Pendry [16], Tong [31], and Van Hove [10]. Perhaps the most widely known is the Barbieri–Van Hove symmetrized automated LEED (SATLEED) package and its derivatives (MSATLEED, SATCLEED, ATLMLEED) [32, 33]. Alternatives include TensErLEED [14], described below, as well as DL_LEED [34], LEED90 [35], CAVLEED [36], CLEED [37] and LEEDFIT [38], though the latter four do not implement the tensor-LEED approximation, while DL_LEED has a partial tensor-LEED implementation [39]. The various packages encompass different capabilities, but none automates the structure-optimization process. Recently, the AQuaLEED package was published to simplify use of the SATLEED codes through a Python wrapper program [40, 41] in a similar, though much less broad, approach as the ViPErLEED package described here. A Python wrapper has also been recently developed for CLEED [42].

The theoretical foundation and code architecture of the TensErLEED package are described in detail in Ref. 14. Briefly, the package consists of three main modules, shown in Fig. 3(a), which perform different segments of the LEED-$I(V)$ calculation. The “reference-calculation” module performs full-dynamic LEED-$I(V)$ calculations, producing the complex amplitudes $A_{\mathbf{g}}$ of the diffracted beams, $\mathbf{g}$, as implemented previously in the Van Hove–Tong code [10]. In addition, it also stores information necessary for the tensor-LEED approximation in separate output files, named “Tensors” in Ref. 14. The second TensErLEED segment is a “delta-amplitudes” calculation, which computes the changes $\delta{}A_{i,\mathbf{g}}$ for a small perturbation to the position, vibration amplitude, or site occupation (including mixing of chemical species as described by the average $t$-matrix approximation [43]) of an atom $i$. These $\delta{}A_{i,\mathbf{g}}$ are calculated for all atoms under variation and for each step in a given atom’s range of displacements. The third segment, the “search”, aims for the set of $\delta{}A_{i,\mathbf{g}}$ that, together with the reference amplitudes $A_{\mathbf{g}}$, minimizes the $R$ factor (i.e., the disagreement) between the calculated and experimental beams. Simply put, the full-dynamic reference calculation is a “precise” LEED calculation for a given reference structure, while the $\delta{}A_{i,\mathbf{g}}$ calculation and “search” compose the structural optimization using tensor LEED.

Altogether, the TensErLEED package [14] is a powerful and versatile tool for LEED-$I(V)$ calculations, but it has shortcomings. Employing the package with its custom settings demands significant experience and poses a serious burden of entry for LEED-$I(V)$ experimenters.

ViPErLEED drastically simplifies the user experience and expands the functionality and performance of existing solutions. The Python package controlling the computational aspects is named viperleed.calc. It manages the LEED-$I(V)$ calculations — acting as a wrapper for the established TensErLEED modules — and contains higher-level logic and automation. By defining suitable defaults based on the theoretical structure and experimental data, viperleed.calc dramatically cuts down the required user input. Moreover, it provides additional functionality, such as optimization procedures not supported in the tensor-LEED formalism, automatic detection and preservation of surface and bulk symmetries during structural optimization, and acceleration of calculations through automated search procedures and parallelization.

viperleed.calc is designed with the vast scientific Python ecosystem in mind and leverages established numerical and scientific libraries such as NumPy [44], SciPy [45], Matplotlib [46], and Scikit-learn [47]. At the same time, it has an open application programming interface for future projects to exploit its novel features, such as the automated plane-group symmetry detection and structure symmetrization introduced in Section II.4.2.

LEED-$I(V)$ users wishing to work with existing solutions like TensErLEED, mostly based on the Van Hove–Tong codes [10], are faced with several challenges. The first hurdle regards the input files: Specific input formats are expected, and these may be difficult for new users to understand. Specifically, standardized structure files, often obtained as a result of density-functional theory (DFT) calculations, cannot be used as input directly. Instead, they have to be translated to “LEED formats” manually, which is tiresome and error prone. The standard output files are also problematic. To the best of the authors’ knowledge, none of the existing codes applying tensor-LEED calculations output the optimized structure in a standard crystal-structure format. As an example, TensErLEED optimization returns only a list of indices corresponding to the displacements of individual atoms that minimize the $R$ factor [14]. Translating these raw indices back to atomic coordinates of the best-fit structure is, again, cumbersome and a likely source of mistakes. When multiple optimization steps are required, all of these steps must be performed by hand repeatedly. Lastly, users are often required to manually compile the source code at runtime. This is for example the case for TensErLEED: many array dimensions depend on the input data and Fortran 77 does not allow dynamic memory allocation.

viperleed.calc addresses these issues. It accepts input in a standard format and automatically translates it into the custom format required by TensErLEED. Likewise, TensErLEED outputs are translated back to the same standard formats by keeping track of the inputs and variations. viperleed.calc also removes the need for manual compilation at runtime. For each TensErLEED segment, it performs the steps shown in Fig. 3(b): (i) Translate the user input into the format required by TensErLEED, (ii) compile the relevant TensErLEED code with appropriate array dimensions, (iii) execute the compiled code as a subprocess with adequate parallelization and with real-time monitoring if required, and (iv) process the TensErLEED output, producing user-readable output or input for the next TensErLEED segment. The compilation and execution of source code is a potential safety threat, as the source files may be modified with malicious intent. To ensure that only the expected TensErLEED code is executed, viperleed.calc contains hardcoded checksums for all expected Fortran files, and by default will only compile and run files conforming to these checksums.

Finally, ViPErLEED also includes a “bookkeeping” utility that keeps track of the performed calculations. It automatically updates a plain-text file (history.info) with information about the run segments, the resulting $R$ factors, and optional notes by the user. Additionally, it archives input and output files for every calculation into a history directory.

The main input parameters defining viperleed.calc behavior are passed via the PARAMETERS input file. In the minimal case, this file defines which segments of TensErLEED should be executed and provides some information on how to interpret the structure files. However, the PARAMETERS file also allows in-depth control of TensErLEED calculation parameters and ViPErLEED defaults. More information on these parameters and the choice of defaults will be given in Section II.4. A comprehensive list of parameters and their default values is found in the ViPErLEED documentation [48].

Other than the general run parameters, the only input strictly required for calculating $I(V)$ curves is the structural data, that is, an atomic model for the diffracting surface. The chosen standard is the POSCAR format used by the Vienna ab-initio simulation package (VASP) [49]. This choice is based on the assumption that many LEED-$I(V)$ calculations start from a DFT-optimized structure; POSCAR files can easily be generated by standard structure viewers such as VESTA [50]. However, structural data can also be passed to ViPErLEED directly from the atomic simulation environment (ASE) package [51], which supports a variety of structural formats. Note that the surface is required to face the $+z$ direction, with the first two unit-cell vectors $\mathbf{a}$ and $\mathbf{b}$ in the surface $(x,y)$ plane. Additionally, viperleed.calc assumes the bottommost layers (lowest $z$ coordinates) to be bulklike. This means that symmetric slabs commonly used in DFT are not allowed. However, a ViPErLEED-compatible slab can be obtained from a standard “DFT slab” by cutting below the “fixed” bulklike layers.

If structural optimization is to be performed, viperleed.calc requires two additional inputs: experimental $I(V)$ curves (file EXPBEAMS.csv), and instructions on optimization parameters and ranges (file DISPLACEMENTS). For $I(V)$ curves, ViPErLEED uses standard CSV formatting for all input and output. The DISPLACEMENTS file defines the range over which each structural parameter is varied during tensor-LEED optimization, and can contain instructions about the order in which the parameters are optimized. Structural parameters that can be varied are (i) atom positions, (ii) vibration amplitudes, and (iii) site occupation (including changes of statistical chemical composition through, e.g., vacancies, dopants, or alloy constituents). For convenience, the DISPLACEMENTS file supports a variety of shorthand forms to address groups of atoms by element, layer, specific site type, or atom number. The exact syntax for the DISPLACEMENTS file is described in the ViPErLEED documentation [48]. Plane-group symmetries of the surface are automatically preserved by viperleed.calc during displacements, such that symmetry-equivalent atoms in the input structure are kept symmetry equivalent during the optimization, unless the user explicitly reduces the symmetry by setting the desired plane group in the PARAMETERS file (see Section II.4.2).

Besides the mandatory input discussed above, users can optionally provide additional files. Otherwise, these are automatically generated by viperleed.calc:

• •

  First, an IVBEAMS file that specifies the $I(V)$ curves to be calculated. This can be automatically generated to match available experimental data. When no experimental data is given, IVBEAMS must be supplied by the user.

• •

  Second, a starting guess for vibration amplitudes [52] and site occupations. This information can be included in the VIBROCC file. Otherwise, sites are assumed to be “chemically pure”. If vibration amplitudes are not specified, an initial guess can be generated based on the Debye temperature $\Theta_{\text{D}}$ of the material and the experimental sample temperature $T$ — both specified in PARAMETERS (see the ViPErLEED documentation [48] for more details). In addition, a scaling factor may be defined for groups of atoms to reflect differences in the coordination environments. For example, undercoordinated atoms at the surface are expected to have a larger vibration amplitude. Since the vibration amplitudes themselves are parameters in the tensor-LEED optimization, this approximation is sufficient to obtain a reasonable starting point even when the Debye temperature is not known precisely.

• •

  Finally, energy-dependent phase shifts for electrons scattering on each type of atom in the structure. If not provided by the user (via a PHASESHIFTS file), these phase shifts will be generated automatically based on the input structure, using the EEASiSSS program [25], distributed with the ViPErLEED package with permission by the author.

ViPErLEED output files use the same formats as the input files. This ensures that, if further optimization steps are required, output from previous ViPErLEED runs can be used directly as input for subsequent runs. The optimized structure is written to the files POSCAR_OUT and VIBROCC_OUT, and calculated $I(V)$ curves are again written in CSV format. In addition, viperleed.calc generates PDF files with plots comparing the calculated beams to experimental data, annotated with $R$ factors for each beam. Some plots as output by viperleed.calc for the $\alpha\text{-Fe}_{2}\text{O}_{3}(1\overline{1}02)$-$(1\hskip 1.00006pt\times\hskip 1.00006pt1)$ example system are shown in Fig. 4. More examples can be found in the ViPErLEED documentation [48]. During structural optimization, a PDF file illustrating the search progress and state of convergence is generated in addition to the text-based log files (see Fig. S1 of the Supplemental Material [27]).

For diagnostic purposes, several supplementary files are written during each viperleed.calc run and stored in a directory called SUPP. These include the inputs and outputs of the TensErLEED programs, as well as modified POSCAR files generated during the structure analysis for troubleshooting.

As discussed in Section II.2, one goal of ViPErLEED is to reduce the number of parameters required as user input. While some TensErLEED parameters can simply be set to a default value that is normally appropriate, others must be chosen in accordance with the structure in question. However, in most cases, it is still possible to define case-specific defaults through educated guesses.

To give an example, the scattering phase shifts $\delta_{\ell}(E)$ depend on the angular momentum $\ell$. At large $\ell$, the phase shifts become small. Since the computing time strongly depends on the $\ell$ values (Section I.3), a cutoff $\ell_{\text{max}}$ must be defined. viperleed.calc automatically selects the cutoff such that phase-shift values at $\ell>\ell_{\mathrm{max}}$ are smaller than a user-defined threshold (parameter PHASESHIFT_EPS). Similar considerations apply to many parameters that can be determined from the structural input. For most low-to-medium-complexity cases, viperleed.calc requires user input of less than ten parameters for a standard optimization. A full list of parameters, their defaults, and descriptions on how they are determined is given in the ViPErLEED documentation [48], and some example inputs are included in the Supplemental Material [27]. Here, the focus is on two non-trivial sets of parameters that viperleed.calc detects from the structural input: the repeat unit of the bulk, and plane-group symmetry.

A typical LEED-$I(V)$ optimization with TensErLEED proceeds as follows: (i) execution of a full-dynamic calculation for the starting guess structure; (ii) optimization along a given displacement direction; (iii) further optimizations along other directions. [Steps (ii) and (iii) are distinct because the TensErLEED search module only supports 1D displacements for atom positions, such that multiple executions are required to fully optimize 3D atomic positions.] Steps (i) and (iii) are repeated until convergence. In TensErLEED, the user must perform all of these steps manually. This becomes especially cumbersome as one segment’s output must be manually translated into an input for the next segment. viperleed.calc simplifies the workflow. First, it handles the input/output processing internally. Second, it allows arbitrary queuing of segments with the RUN parameter. Finally, it allows to specify multiple sets of displacements in the DISPLACEMENTS file, such that they are automatically executed in the requested order. These optimizations can be looped until the $R$ factor converges, allowing for an autonomous search in a wide parameter space. The ViPErLEED documentation [48] gives more details on the RUN parameter and DISPLACEMENTS file. As described in Section S1 of the Supplemental Material [27], viperleed.calc also speeds up each optimization run in TensErLEED by monitoring the search progress and adjusting convergence parameters accordingly.

Besides automating the optimization workflow, viperleed.calc also enhances the parallelization of LEED-$I(V)$ calculations. The parallelization offered by the native TensErLEED code is limited to the “search” module. viperleed.calc introduces additional parallelization at the Python level. Parallelization is applied to the full-dynamic calculation by assigning individual energies to separate processes. This is possible since each energy step can be calculated independently from the others. This approach enables further acceleration through a dynamic use of $\ell_{\text{max}}$ at each energy: At low energies, where phase shifts of large angular momenta are negligible, small $\ell_{\text{max}}$ values are used. Parallelization is also applied to the $\delta{}A_{i,\mathbf{g}}$ calculation. Here, the atoms under variation are independent and can be calculated in a parallel manner. Note that, since the single calculations are vastly more computationally expensive than the overhead, the parallelization offered by viperleed.calc comes close to the efficiency theoretically achievable with parallelization directly in the Fortran code.

Although ViPErLEED is backwards compatible with the most recent version of TensErLEED at the time of writing (version 1.6), it was expedient during development to also make some improvements on the TensErLEED codebase. This resulted in a new version 2.0 of TensErLEED [57], which is released with ViPErLEED but can also be used independently. Changes include some minor bug fixes, increased precision of some computational variables, and improved computational performance, especially of the full-dynamic “reference” calculation. The speed-ups mainly come from replacing custom matrix-algebra routines by better-optimized BLAS and LAPACK versions [58], skipping redundant calculations, and refactoring certain routines for more efficient memory access.

The user interface (i.e., format of the input and output files) in TensErLEED 2.0 was modified to extend the size limits of systems that can be calculated. Previous TensErLEED versions severely limited the size of the unit cell, as well as related quantities such as the number of diffracted beams to calculate. The TensErLEED 2.0 input now allows larger values. Furthermore, some variables that were obsolete in the TensErLEED code but were still required in the input were removed, and others (e.g., imaginary part of the inner potential VPI and the filament work function WORKFN) were moved from hardcoded values in a runtime-created Fortran script (muftin.f) to standard input parameters.

Modifications to the TensErLEED code by the Solid-State Physics group in Erlangen and the ViPErLEED team are released together with ViPErLEED under the terms of the GNU General Public License version 3 (or any later version) [59].

The ViPErLEED version 1.0 released with this work is a fully functional open-source package, licensed under the GNU General Public License version 3 or later [59]. Further development is in progress and will be published on the GitHub repository [60].

Already in version 1.0, ViPErLEED can start calculations based on structure input in the form of ASE Atoms objects [51]. This interface will be expanded to allow full control of all ViPErLEED features. This should, for the first time, allow easy high-throughput LEED-$I(V)$ calculations for automatically generated surface structures. Many computational groups in the field are currently striving to solve complex surface structures without relying on human intuition [61, 62, 63, 64, 65, 66, 67]. As demonstrated in Ref. 68, the $R$ factor can be used as a cost function for guiding machine-learned structure prediction. The sensitive feedback provided by LEED $I(V)$ may provide an attractive and new use of the technique for structure search in conjunction with total energies yielded by DFT calculations. A major benefit of using LEED $I(V)$ for structure determination is its direct link to experimental data.

Another major change currently under development concerns a replacement for the TensErLEED optimization (“search”) algorithm. While ViPErLEED will continue to use the TensErLEED code for the full-dynamic calculations and calculations of amplitude changes $\delta{}A_{i,\mathbf{g}}$, the TensErLEED algorithm for finding the combination of $\delta{}A_{i,\mathbf{g}}$ that minimizes the $R$ factor is somewhat dated. In the future, the current one-dimensional, grid-bound search may be replaced by a continuous search by calculating $\delta{}A_{i,\mathbf{g}}$ values on demand. This would allow direct application of, and free choice among established global optimization algorithms available from many optimized Python packages (e.g., SciPy [45], Scikit-Optimze [69], Keras [70], Ray [71], JAX [72]). Most of these modern evolutionary or (stochastic) gradient-based algorithms are likely to be much more effective at optimizing the structure than the original custom algorithm implemented in TensErLEED.