Hydrogen atom scattering at the Al₂O₃(0001) surface: a combined experimental and theoretical study

Abstract

Investigating atom–surface interactions is the key to an in-depth understanding of chemical processes at interfaces, which are of central importance in many fields – from heterogeneous catalysis to corrosion. In this work, we present a joint experimental and theoretical effort to gain insights into the atomistic details of hydrogen atom scattering at the α-Al₂O₃(0001) surface. Surprisingly, this system has been hardly studied to date, although hydrogen atoms as well as α-Al₂O₃ are omnipresent in catalysis as reactive species and support oxide, respectively. We address this system by performing hydrogen atom beam scattering experiments and molecular dynamics (MD) simulations based on a high-dimensional machine learning potential trained to density functional theory data. Using this combination of methods we are able to probe the properties of the multidimensional potential energy surface governing the scattering process. Specifically, we compare the angular distribution and the kinetic energy loss of the scattered atoms obtained in experiment with a large number of MD trajectories, which, moreover, allow to identify the underlying impact sites at the surface.

---

1 Introduction

Heterogeneously catalyzed reactions are of utmost economical importance, as a substantial part of the bulk products in chemical industry is made through reactions at solid catalyst surfaces.^1,2 Prominent examples are the Haber–Bosch process for ammonia synthesis,³ the cracking of long-chain hydrocarbons in the petroleum industry,⁴ and the oxidation of sulphur dioxide to sulphur trioxide finally yielding sulphuric acid.^5,6 The improvement of such catalytic systems crucially depends on a detailed understanding of the reaction mechanisms at the atomic scale that are governed by the potential energy surface (PES), which is a high-dimensional function providing the energy and forces for any given atomic configuration. Thus, the PES determines the dynamics of the reactants at the surface, the making and breaking of bonds and, finally, product formation.

Insights into the interactions between atoms or molecules and surfaces can be obtained from both, experimental as well as theoretical studies. In particular atom and molecular beam experiments under well-defined conditions can provide invaluable details of processes at single-crystal surfaces.^7–9 For instance, information about the interaction potential and energy transfer between atoms or molecules and surfaces can be derived from inelastic scattering experiments. Moreover, using rare-gas atoms, surface phonons of a material can be studied.^10,11 Further, the role of rotational, vibrational and translational energy on surface reactivity is of central interest in molecule–surface scattering experiments.^8,12–15 Finally, scattering of open–shell atoms like hydrogen or oxygen provides detailed information about chemical bond formation in surface reactions.^9,16

Next to experiment, theoretical studies have become an indispensable tool in surface science and heterogeneous catalysis in the past two decades.^17–20 Still, in spite of the substantial increase in computational power, first-principles methods like density functional theory (DFT), which have been very successfully applied to the calculation of properties like adsorption energies and surface structures, remain too demanding for determining a large number of molecular dynamics (MD) trajectories when computing the forces on-the-fly for each new MD step. Thus, it is common practice to perform extended MD simulations for scattering at surfaces employing atomistic potentials, which – after fitting to accurate electronic structure data – can provide a direct relation between atomic positions and forces while they are much cheaper to evaluate. A wide range of methods has been proposed to represent atomistic potentials for surface scattering processes,^21–23 but constructing high-dimensional PESs of first-principles quality taking all the surface degrees of freedom explicitly into account has remained a substantial challenge.

In recent years, the introduction of machine learning potentials (MLP)^24–29 has led to a paradigm change in the development of PESs, since flexible machine learning algorithms allow to combine the accuracy of electronic structure methods with the efficiency of simple empirical potentials. In fact, molecule–surface interactions have been among the early applications of MLPs,^30–34 which in the first years still employed a frozen-surface approximation to restrict the complexity of the PES. Starting with the introduction of high-dimensional neural network potentials (HDNNP) in 2007,^35,36 the construction of MLPs for high-dimensional condensed systems has become possible. Nowadays, many different types of MLPs are available,^37–42 which allow to explicitly take all degrees of freedom into account even for systems containing thousands of atoms. This methodical progress now enables to study atom- and molecule–surface interactions with full dimensionality including a mobile surface, and numerous applications have been reported to date.^43–47

In spite of these advances, so far most studies have addressed the surfaces of rather simple materials like metals,^44,48–50 solid xenon⁵¹ or graphene,^43,52i.e., pure elemental systems, while scattering at binary compounds like oxides or thin surface oxides⁵³ has been rarely studied to date. Only recently, a first H-atom beam experiment has been reported for the interaction with α-Al₂O₃.⁵⁴ Its surface termination is well-studied under high vacuum conditions and does not show any surface reconstruction. The material is also often applied, either as catalyst⁵⁵ or as supporting material,⁵⁶ increasing the need for a deeper understanding of processes at its interfaces.

Here, we report a combined theoretical and experimental study addressing the interaction of H-atoms with the α-Al₂O₃(0001) surface with the aim to unravel the atomistic details of the scattering process. For this purpose we perform H-atom beam experiments employing four different incident conditions to probe the scattering angle distribution as well as the kinetic energy loss of the atoms leaving the surface. Unlike other systems, such as those involving vibrations or steric effects of scattering molecules, in the present work on an atomic beam only the translational energy is involved. We study the effect of this energy for two different incidence angles. For a deeper analysis, in parallel, we carry out large-scale molecular dynamics simulations based on energies and forces obtained from a HDNNP, which after training provides a first-principles quality description of the scattering process. This combination of methods is used for a detailed analysis of the quality of the PES for different initial beam kinetic energies and incidence angles, which allows us to assess the reliability of the employed exchange correlation functional for scattering at different α-Al₂O₃ surface sites.

2 Methods

2.1 Experimental methods

The hydrogen atom scattering instrument employed in this work is described in detail in ref. 9 and 57. It utilizes a mono-energetic H-atom beam formed by photolysis of a supersonic jet of hydrogen iodide molecules using an excimer laser operating with KrF (248 nm). The generated H-atoms pass through a skimmer and two apertures separating two differential pumping regions and hit the sample in the ultrahigh vacuum (UHV) scattering chamber. The sample is mounted on a 6-axis manipulator that allows the variation of the polar and azimuthal incidence angles. Scattered H-atoms are detected by Rydberg atom tagging time-of-flight (TOF),⁵⁸ where two laser pulses excite the H-atoms to a long-lived (n = 34) Rydberg state. The neutral Rydberg atoms travel 250 mm, pass a grounded mesh and are field-ionized and detected with a multichannel plate (MCP) detector. A multichannel scaler records the TOF distributions that are converted to kinetic energy distributions. The detector is rotatable in the plane containing the H-atom beam and the surface normal, providing TOF spectra at various scattering angles. The alumina sample was cleaned by annealing for several hours in an oxygen atmosphere (10 × 10⁻⁶ mbar) at 600 °C. The cleanliness and structure of the surfaces are monitored by Auger electron spectroscopy (AES) and low-energy electron diffraction (LEED). Based on AES and LEED we conclude that the surface has a (1 × 1) structure with Al termination, Al–O₃–Al trilayers^59–62 and a step density of about 2–4%. Moreover, it is important to note that depending on the stacking of the trilayers there are two possible terminations related by a mirror operation resulting in two possible incident azimuthal angles ϕ_i with respect to the [101̄0] surface direction of either 0° or 180°⁶³ that cannot be distinguished with the methods available in experiment, which needs to be considered in the simulations. Further details about the sample characteristics are given in the ESI.† The inelastic scattering of the H-atoms was performed for initial kinetic energies of 0.99 eV and 1.92 eV and incidence polar angles of θ_i = 40° and θ_i = 55° resulting in four different sets of scattering conditions.

2.2 High-dimensional neural network potentials

The full-dimensional potential energy surface of the hydrogen atom scattering process at α-Al₂O₃ is represented by a high-dimensional neural network potential. HDNNPs are a frequently used type of machine learning potential introduced by Behler and Parrinello in 2007.³⁵ In second-generation HDNNPs as employed here,⁶⁴ the total energy of the system is constructed as a sum of atomic energies,that depend on the local atomic environments up to a cutoff radius R_c, which has to be chosen large enough to include all relevant atomic interactions. The positions of all neighboring atoms inside the resulting cutoff spheres are described by vectors of atom-centered symmetry functions (ACSF),⁶⁵ which are invariant with respect to rotation, translation and permutation of the system and in combination with eqn (1) allow to explicitly include all atomic degrees of freedom. The atomic energies are then computed as outputs of element-specific atomic neural networks as a function of the respective input ACSF vectors. Since only the Cartesian coordinates of all atoms in the system are required for the calculation of the ACSFs, HDNNPs are able to describe the making and breaking of bonds, which is an essential requirement for studying scattering at surfaces.

The weight parameters of the atomic neural networks are obtained in an iterative training process making use of energies and forces obtained from electronic structure calculations.^66,67 Since atomic energies are not quantum mechanical observables, total energies are used as target properties, and the partitioning into atomic energies is done automatically during the training process. The forces, which are needed for training as well as for running molecular dynamics simulations, can be computed as analytic derivatives of the HDNNP energy. More details about HDNNPs, the underlying methodology and typical applications can be found in several reviews.^36,64,66,68

3 Computational details

3.1 Density functional theory calculations

The DFT reference energies and forces for training the HDNNP have been calculated using the Fritz-Haber-Institute ab initio molecular simulations (FHI-aims) code (version 171221_1)⁶⁹ applying the “light” settings for the integration grid and for the basis sets, which consist of numerical atomic orbitals. Spin-polarized calculations were performed for structures containing hydrogen atoms, while α-Al₂O₃ bulk and slab structures without hydrogen atoms have been treated as closed–shell systems in a numerically consistent way. Γ-Centred k-point grids have been used for all periodic systems. For the bulk α-Al₂O₃ structures containing 40 atoms a 6 × 6 × 2 k-point grid was used, while for the (2 × 2) slab supercells a 2 × 2 × 1 k-point grid has been employed. Using these k-point grids, the formation energy of bulk α-Al₂O₃ and the adsorption energies of hydrogen atoms at different surface sites are converged to about 0.5 meV. The convergence criteria for the electronic self-consistency of the single point calculations have been set to 10⁻⁶ eV for total potential energies and 10⁻⁴ eV Å⁻¹ for the forces.

Several exchange correlation functionals have been tested in combination with Tkatchenko–Scheffler dispersion corrections,⁷⁰ including the GGA functionals PBE⁷¹ and RPBE⁷² as well as the PBE0⁷³ hybrid functional. Fig. S3 in the ESI† shows two one-dimensional cuts through the potential energy surface for a hydrogen atom approaching the surface on top of an aluminium and an oxygen surface site, respectively. For both sites the observed interaction energies at the potential minimum are within a range of 0.1 eV for all functionals, which is the expected uncertainty with respect to the description of exchange and correlation. This is about one order of magnitude smaller than the lower initial kinetic energy of 0.99 eV used in experiment, even if a potential energy amount of 0.1 eV would be converted entirely to H-atom kinetic energy. At much lower experimental kinetic energies, however, the accuracy of the employed functional would become highly relevant. At both surface sites the PBE functional is more attractive than RPBE with the minimum region showing no clear trend for a better agreement with the hybrid functional. Hence, in the present work, which strongly depends on the description of the repulsive walls due to the high kinetic energies of the H-atoms, we will employ the RPBE functional that is in reasonable agreement with PBE0 for smaller, i.e., repulsive, atom–surface distances. Performing a very large number of PBE0 calculations would be very demanding for the construction of the training set while on the other hand using this hybrid functional does not yield a qualitatively different description of the PES.

3.2 Construction of the reference set

The reference data set contains bulk geometries and slab structures of the α-Al₂O₃(0001) surface, both with and without an H-atom at various positions above the slab. The bulk structures contain one unit cell of α-Al₂O₃, while for the slabs (2 × 2) supercells of the α-Al₂O₃(0001) surface unit cell have been constructed with the vacuum aligned in z direction to avoid lateral interactions between the periodic images of the H-atoms. The slabs consist of 8 Al–O₃–Al trilayers as shown in Fig. 1. The surfaces are terminated by one Al-atom per unit cell on top of an oxygen layer, which is called the “half metal” termination and corresponds to the surface structure present in experiment.^60–62 The atoms of the bottom four trilayers have been frozen at their bulk equilibrium positions, while the atoms in the top four layers have been relaxed and are also mobile in the MD simulations, which ensures that the frozen surface atoms in the deep layers are outside the cutoff spheres of the ACSFs describing the hydrogen atom environment. Utilizing an analysis of the atomic interactions based on the Hessian matrix of the slab,⁷⁵ we found that the remaining interactions between the deep frozen layers and the H-atom are negligibly small even for H-atom positions close to the surface (see ESI†). The vacuum between the slabs was set to 21 Å. At a distance from the surface of about 6.3 Å corresponding to the ACSF cutoff only very small forces of roughly 0.003 eV Å⁻¹ act on the hydrogen atoms such that atomic interactions at larger separations can be neglected to a good approximation. The functional form of the cutoff function⁶⁵ ensures a smooth decay of the forces to zero at the cutoff radius in the HDNNP. The size of the vacuum prevents any notable interactions between the H-atom and the opposite side of the slab in the periodic DFT calculations and the HDNNP even for the largest H-atom-surface distance of 6.8 Å we studied.

Fig. 1 Side (a) and top view (b) of the relaxed (2 × 2) α-Al₂O₃(0001) supercell slab. Al atoms are shown in grey, O atoms are displayed in red. The white atoms represent the regular H-atom grid employed for the initial screening of the potential-energy surface using a frozen surface. The orange lines in (b) represent the boundaries of the symmetry-unique wedges, and two examples for the three-fold rotational symmetry axes aligned perpendicular to the surface are highlighted. The figures have been created using Ovito version 3.8.4.⁷⁴

The reference data points were generated combining various approaches. Initially, 133 bulk and 200 clean slab structures were generated with random atomic displacements of up to 0.2 Å. Initial HDNNPs based on these structures were used in MD simulations in the canonical ensemble at increasing temperatures up to 400 K to identify important missing surface structures through active learning.^76–78 Structures showing an energy or force variance larger than 2.8 meV atom⁻¹ or 0.5 eV a₀⁻¹ (with a₀ being the Bohr radius) for different preliminary HDNNPs, respectively, were added to the data set.

Hydrogen atom positions above the surface were initially included using an equidistant grid with spacing of 0.5 Å in the x and y Cartesian directions, starting at the vertical position of the topmost oxygen atom layer and up to 6.75 Å above this layer within the symmetry-unique part of the unit cell (Fig. 1) employing the fixed structure of the clean relaxed surface. In the z direction the distances between the H-atom positions are increasing with distance from the surfaces due to the reduced variance in the atomic interactions at larger distances. Up to a distance of 1.0 Å from the topmost oxygen layer a vertical grid spacing of 0.125 Å has been employed, followed by a vertical spacing of 0.25 Å for distances up to the limit of 6.75 Å. From this grid, all strongly repulsive structures with a distance between the H-atom and any surface atom below 0.5 Å have been discarded. In the next step, active learning has been performed employing molecular dynamics simulations of the scattering process using the experimental incidence conditions including a mobile surface to sample the surface degrees of freedom. Additional structures obtained in these trajectories were additionally selected by farthest point sampling⁷⁹ as discussed in the ESI† to improve the representation of rarely visited geometries.

3.3 Construction of the high-dimensional neural network potential

For the construction of the HDNNP the RuNNer code (version from August 22, 2019) was used.^36,66 The atomic neural networks consist of two hidden layers containing 19 neurons each with n^α_G input neurons for atoms of element α corresponding to the respective number of ACSFs, and one output node providing the atomic energy. The ACSF cutoff radius is R_c = 12a₀ (6.35 Å), which is sufficient to include a major part of the atomic interactions as discussed in Section 3.2. A list of the employed atom-centered symmetry functions can be found in the ESI.† Since there is only one H-atom per structure, ACSFs for the description of H–H interactions have not been included. For the construction of the HDNNP the atomic energies were removed from the total energies before training such that numerically more favorable binding energies are learned. The HDNNP was then trained to reproduce these binding energies and the DFT atomic force components of the reference structures. The reference data set was randomly split into a training set consisting of about 90% of the structures used for adjusting the weight parameters and a testing set of the remaining 10% of the structures for assessing the quality of the HDNNP for unknown structures. The parameters of the Kalman filter⁸⁰ were set to λ = 0.98000 and ν = 0.99870. Overall, the structures cover an energy range of 0.446 eV atom⁻¹. Atoms with forces up to 13.6 eV Å⁻¹ have been included in the data set for monitoring the training errors, but only forces up to 3 eV Å⁻¹ have been used to update the weights in the training process.

3.4 Molecular dynamics simulations

To run MD simulations the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) (version from 16th March 2018)⁸¹ including the n2p2 extension for HDNNPs⁸² was used. MD simulations for active learning were run in the canonical NVT (bulk and clean surfaces) and in the microcanonical NVE (trajectories for H-atom scattering) ensemble with a time step of δt = 0.5 fs for simulations without a hydrogen atom and δt = 0.25 fs for simulations including a hydrogen atom. The velocity Verlet algorithm was used as integrator.⁸³ To simulate the NVT ensemble a Nosé–Hoover thermostat⁸⁴ with a temperature damping parameter of 0.01 ps was used.

For the MD trajectories of H-atom scattering, i.e., slab structures at room temperature, atomic velocities of the surface atoms were taken from NVT trajectories at 300 K in intervals of 0.1 ps. Then, the H-atom was placed at a random lateral position 6.3 Å above the surface with velocity vectors corresponding to the experimental kinetic energy, polar and azimuthal angle. The trajectories were run in the NVE ensemble and terminated after 0.5 ps or before, when the H-atom surface distance exceeded a value of 7.8 Å. Trajectories were considered in the analysis of the kinetic energy loss and the angular distributions if the H-atom was at least at a distance of 7.3 Å from the slab at the end of the trajectory. The resulting kinetic energies and scattering angles of the H-atoms were calculated using the velocity vector of the last time step of the trajectory. The detector counting the H-atom flux is defined as a vector with the specified polar angle θ_s towards the surface normal and in plane with the incident azimuthal angle ϕ_i to mimic the position of the detector in the experiment. The experimental detection angle γ, which is the angle between the velocity vector of the H-atom and the detector, is about 1.5°, while in the MD simulations larger angles up to 5° have been included to improve statistics, which has very little effect on the obtained results (for a discussion see ESI†).

3.5 Phonon calculations

Phonon band structures were calculated with the phonopy code (version 2.17.1)^85,86 using both, force constants from FHI-aims RPBE DFT calculations as well as force constants calculated using the HDNNP. The force constants were determined using a geometry optimized primitive α-Al₂O₃ cell containing 10 atoms. The geometry optimization using the HDNNP was carried out using LAMMPS⁸¹ with the n2p2⁸² extension. To determine the correct symmetry, the symmetry precision threshold was set to 5 × 10⁻⁴. The interface between phonopy and RuNNer was created using a Python script.

4 Results and discussion

4.1 Validation of the HDNNP

The final data set obtained after several active learning cycles consists of 15 812 structures, which include 808 bulk configurations and 15 004 slab geometries (12 704 with H-atom and 2300 clean surface structures). The final HDNNP has a root mean square error (RMSE) of 0.746 meV atom⁻¹ for the energy and 0.103 eV Å⁻¹ for the atomic force components for the test data set, which is in the typical order of magnitude of state-of-the-art MLPs. For the training set the RMSE values are 0.257 meV atom⁻¹ and 0.111 eV Å⁻¹, respectively. Therefore, the HDNNP errors with respect to the reference method are significantly smaller than the uncertainty with respect to the choice of the exchange correlation functional. More detailed information about the energy and force errors can be found in the ESI.†

As monitoring the RMSEs of the energies and forces alone is insufficient to judge on the quality of a potential, in a first validation step we calculated the lattice constants and bulk modulus B₀ for bulk α-Al₂O₃. As can be seen in Table 1, the HDNNP predictions are in excellent agreement with the underlying RPBE-DFT data, while compared to experiment the RPBE functional slightly overestimates the lattice constants and underestimates the bulk modulus. In a next step, we have computed the phonon band structure of bulk α-Al₂O₃ using the phonopy program.^85,86 The phonon spectrum obtained from the HDNNP in Fig. 2 shows excellent agreement with DFT in particular at low frequencies. Overall, the phonon band structure, which is important for energy dissipation after the scattering process, is well represented, and we conclude that the HDNNP provides a reliable description of bulk α-Al₂O₃.

Table 1 Calculated and experimental lattice constants a and c and bulk modulus B₀ of bulk α-Al₂O₃⁸⁷

Parameter	RPBE	HDNNP	Exp.
a/Å	4.827	4.824	4.7554^87
c/Å	13.145	13.139	12.991^87
B 0/GPa	222	216	253^88

---

Fig. 2 Phonon band structure of bulk α-Al₂O₃ computed for the HDNNP and DFT employing the RPBE functional using the phonopy program.^85,86

In order to ascertain the quality of the structural description of the clean surface, we compared the relaxed interlayer distances of the clean surface optimized with RPBE DFT and with the HDNNP. When optimizing the surface with DFT, the topmost three interlayer distances relax from the initial bulk separations of 1.59 Å (Al–O₃), 1.59 Å (O₃–Al) and 0.96 Å (Al–Al), respectively, to 0.18 Å (Al–O₃), 1.69 Å (O₃–Al) and 0.52 Å (Al–Al). Consequently, there is a substantial inwards relaxation in particular of the top Al layer, and in the HDNNP relaxation this Al–O₃ layer distance is about 0.15 Å with an error of only about 0.016 Å. For all other interlayer distances the differences between the HDNNP and RPBE results are below 0.01 Å.

As an initial assessment of the quality of the hydrogen atom–surface interactions, we have analyzed the 3D-PES of a H-atom at different positions above the surface, which, as described in Section 3.2, also formed the starting point of the reference data generation beyond this grid. Fig. 3 shows several 2D cuts through this PES at different H distances from the surface, which is frozen in its relaxed geometry for this purpose. Overall, the RPBE and the HDNNP PESs displayed in Fig. 3(a) and (b) are very similar, as can also be seen in the energy difference plot shown in panel (c). The total energy differences are found to be typically between −0.01 and 0.01 eV, which corresponds to an error of 0.06 meV atom⁻¹. The largest difference is about −0.04 eV in the cuts at 2.0 Å and 2.25 Å, i.e., 0.25 meV atom⁻¹, which is still well below the RMSE of the HDNNP. We conclude that the differences are very small. However, the 3D-PES only probes the interaction between the H-atom and the frozen surface, which is insufficient to check the reliability in MD simulations of the scattering process with fully mobile surface atoms.

Fig. 3 2D cuts at different distances through the PES of the H-atom above a frozen α-Al₂O₃(0001) surface. Panel (a) depicts the (2 × 2) supercell of the RPBE PES while (b) shows the HDNNP PES. Given are relative potential energies ΔE of the entire system with respect to the energy of a H-atom at infinite separation from the surface. The underlying data points correspond to the regular grid of H-atom positions shown in Fig. 1. Panel (c) shows the difference ΔΔE between RPBE and the HDNNP. Note that in contrast to the energy RMSEs these energies are not normalized per atom in the system.

To validate the performance of the HDNNP for full-dimensional trajectories, ab initio MD simulations have been carried out and compared to the HDNNP. The starting structures for the ab inito MD simulations and the atomic velocities for the slab atoms were taken from NVT simulations at 300 K of the clean surface as determined using the HDNNP. The magnitude and direction of the velocities of the H-atoms correspond to the experimental conditions, and we have chosen such conditions that are estimated to result in experimentally detectable scattering events based on preliminary HDNNP trajectories. In total, ten different ab initio MD trajectories have been computed, and the results are shown in Table 2. Since even smallest numerical differences result in diverging MD trajectories, a direct comparison of independently computed ab initio and HDNNP trajectories using the same initial conditions is not possible. Thus, for assessing the reliability of the HDNNP, we have recomputed the structures visited in the ab initio trajectories by single point HDNNP calculations and computed the RMSEs of the energies and the z-components of the forces acting on the H-atom (Table 2). Overall, we find a very good agreement between DFT and the HDNNP for all trajectories with deviations similar to or below the RMSEs of the training and test data sets.

Table 2 Initial conditions of the ab initio MD simulations used to validate the HDNNP, including the incident polar angle θ_i, the incident azimuth angle ϕ_i and incident kinetic energy E_kin,i, as well as the polar angle θ_s and the kinetic energy of the scattered atoms E_kin,s. Moreover, the RMSE values of the total energy and of the z-component of the H-atom force vector have been computed by comparing the ab initio MD trajectory with subsequent HDNNP single point calculations along the trajectories to quantify the deviations between DFT and the HDNNP. The surface has been equilibrated at 300 K

Name	θ i/°	ϕ i/°	E kin,i/eV	θ s/°	E kin,s/eV	RMSE E/meV atom^−1	RMSE FH,z/eV Å^−1
MD1	40	0	0.99	54	0.56	0.39	0.039
MD2	40	0	1.92	25	1.65	0.35	0.041
MD3	40	0	1.92	44	1.69	0.32	0.037
MD4	40	0	1.92	50	1.78	0.38	0.045
MD5	40	180	1.92	−10	1.00	0.34	0.125
MD6	40	180	1.92	32	1.00	0.33	0.047
MD7	40	180	1.92	61	0.83	0.35	0.102
MD8	40	180	1.92	13	0.57	0.47	0.108
MD9	55	0	0.99	49	0.83	0.36	0.022
MD10	55	180	1.92	48	1.70	0.37	0.047

---

Fig. 4 shows a detailed analysis of the scattering process for the example of trajectory MD4 with an incident kinetic energy of 1.92 eV and an incident polar angle of 40°. The trajectory was calculated using RPBE DFT and afterwards the structures were recalculated using the HDNNP. A comparison of the DFT and HDNNP potential energies of the system along the trajectory is shown in Fig. 4(a), while Fig. 4(b) displays a comparison of the DFT and HDNNP z-components of the forces acting on the H-atom. Both, energies and forces, are very well represented by the HDNNP along the entire trajectory, which is shown in a top and a side view in panels (c) and (d). We note that the energy deviations are small in comparison to the initial kinetic energy of the collision of 1.92 eV and thus do not significantly affect the scattering process. A similar agreement has been found also for the other trajectories in Table 2.

Fig. 4 Analysis of ab initio MD trajectory MD4 (see Table 2). Panel (a) shows the DFT and HDNNP potential energies of the structures ΔE_slab+H with respect to the relaxed slab and the H-atom at infinite distance from the surface vs. the time t of the trajectory. The trajectory was calculated using RPBE DFT and afterwards the structures were recalculated using the HDNNP showing very close agreement. Panel (b) shows the z-component of the force vector acting on the H-atom. Panels (c) and (d) show the path (blue) of the H-atom (white) along the surface (top and side views). The surface structure corresponds to the first frame of the simulation, the highlighted H-atom denotes the starting position above the surface.

4.2 Comparison of experimental scattering and simulations

Having validated the accurate description of the DFT-PES by the HDNNP, we now perform large-scale MD simulations and compare the scattering properties of the α-Al₂O₃(0001) surface with experiment using all four scattering conditions as described in Section 2.1. Due to the close agreement between DFT and the HDNNP, a comparison between experiment and simulation will allow to assess the performance of the employed DFT description of the scattering process.

Concerning the experimental side, in view of the rather high complexity of the α-Al₂O₃(0001) surface in comparison to previous atomic beam studies at surfaces of pure elements, we took great care when characterizing the surface prior to the scattering experiments in detail (see ESI†). Still, we cannot exclude that differences between experiment and simulation may be related to experimental shortcomings, and a comparison of experimental and theoretical results needs to be done with care.

In the experiment, incidence polar angle, incidence kinetic energy and surface temperature are well-controlled. Also the incident azimuthal angle is experimentally defined, but here we have to consider that the surface can be terminated in two ways as discussed in Section 2.1. However, both represent very different surface morphologies for a tilted incidence angle of the beam and thus produce significantly different scattering spectra. Due to the lack of more detailed information, for the comparison between experiment and MD we have computed trajectories of both terminations assuming a 1 : 1 ratio in experiment, which in fact also provides the best overall agreement. A further discussion of this aspect can be found in the ESI.†

Fig. 5 shows the experimentally obtained kinetic energy loss spectra for the incidence polar angle of 40° and both initial kinetic energies E_kin = 0.99 eV; 1.92 eV along with the respective data obtained in the MD simulations. The corresponding figure for the incidence polar angle of 55° is given in the ESI.† The distributions have been normalized to a maximum peak height of one. For all incidence conditions a sticking probability of at most 5% has been found (Table 3), which therefore does not significantly affect our results. Overall, it can be seen that the simulations have a slight tendency to overestimate the energy loss of the scattered H-atoms with respect to experiment, achieving better agreement for the larger incidence energy (Table 4).

Fig. 5 Kinetic energy loss distributions for H-atom scattering at (a) E_kin,i = 0.99 eV and θ_i = 40°, (b) E_kin,i = 1.92 eV and θ_i = 40°. The black line shows the experimentally measured distributions, while the red line shows the theoretical distributions. The detected scattering polar angle θ_s is given in each panel. All distributions have been normalized to a maximum of one and the data has been averaged for both azimuthal incident angles of 0° and 180° in the MD simulations.

Table 3 Sticking probabilities P_sticking obtained in the MD simulations for different incidence conditions

E kin,i/eV	θ i/°	P sticking (%)
0.99	40	5.0
0.99	55	5.8
1.92	40	4.1
1.92	55	1.3

---

Fig. 6 compares energy integrated angular distributions of the scattered H-atoms obtained from experiment and simulation. The distributions have been normalized to a maximum of one and the theoretical distributions are again combinations of the two possible surface terminations ϕ_i = 0° and 180°. For all incident conditions the position of the maximum compares well to the experiment and there is very good agreement between experiment and theory for large scattering angles. However, at smaller scattering angles experiment and theory exhibit some deviations, as in experiment scattering angles close to the surface normal are more frequently found compared to the MD simulations.

Fig. 6 Experimental (black) and computed (red) angular distributions of the scattered H-atoms at different incidence polar angles θ_s and kinetic energies. The maxima of the distributions have been scaled to one. The data has been averaged for both azimuthal incident angles of 0° and 180° in the simulations.

4.3 Discussion

Although overall a reasonable agreement has been obtained, there are some deviations between the experimental data and their counterparts determined by molecular dynamics as in the simulations there is a larger kinetic energy loss (Fig. 5) in particular for the lower incidence kinetic energy as well as a narrower angular distribution of the scattered atoms (Fig. 6) compared to the atomic beam results. However, in previous work, good agreement between experiment and simulation was achieved for well-defined, single-element surfaces, i.e., solid Xe,⁵¹ graphene⁴³ and several fcc metal surfaces.⁴⁹ The more pronounced discrepancies in the present study may have several reasons.

First, it has to be ensured that the surface structure used in the simulations is a good representation of the surface that is present in experiment. In case of the α-Al₂O₃(0001) surface, controlling the surface structure is more challenging, as, e.g., there are two possible terminations, which are geometrically inequivalent for the scattering atoms and different ratios of both terminations is essentially unknown. The effect of this ratio on the simulation outcome is discussed in the ESI.† Moreover, the effect of steps and other surface imperfections is not included in the simulations.

Concerning the simulations, the accuracy of the reference electronic structure method poses a limit for the quality of the HDNNP that can be achieved. In the present work a GGA functional has been employed as our initial investigations have shown that for the investigated geometries there are no fundamental differences between the quality of a GGA description and the PBE0 functional in view of the high atomic kinetic energies. However, due to the substantial computational costs of the hybrid functional it has not been possible to construct a full-dimensional HDNNP based on PBE0, and consequently parts of the PES exhibiting larger deviations might have been missed.

We will start discussing the results of the simulations before we address the deviations between experiment and theory in more detail. A first step to gain some basic information about the scattering mechanism is to compare the average energy loss with the predictions of the simple Baule model. The Baule limit describes the maximum energy loss that can be expected in a single-atom surface collision. It is based on a zero impact parameter collision of the impinging H-atom with a resting surface atom. An energy loss significantly larger than the Baule limit is a clear indication for a complex scattering mechanism, for example multiple collisions or energy transfer to several atoms in a single collision. Table 4 compares the Baule limit for the oxygen and the aluminum site to the average energy losses observed in experiment and theory. It can clearly be seen that the simulation, especially for an incidence kinetic energy of 0.99 eV, yields a larger relative energy loss than the simple Baule limit, while in experiment a loss below the Baule limit is observed. This observation indicates that there might be a complex scattering mechanism in the simulations that is not necessarily expected based on the experimental data.

Table 4 Relative average energy losses obtained in MD simulations (sim.) and experiment (exp.) compared to the Baule limit (BL) for θ_i = 40° and θ_s = 50°

E kin,i/eV	Sim.	Exp.	BL, O-site	BL, Al-site
0.99	0.26	0.12	0.22	0.11
1.92	0.22	0.14	0.22	0.11

---

To understand the reason for the large energy loss observed in simulation we took a closer look at the scattering trajectories leading to high energy losses and analysed them in detail. First, we extracted the geometry of closest approach of the H-atom to the surface for each trajectory to allocate the surface impact sites and to correlated them with the energy loss and scattering angle. In Fig. 7 the closest point of approach for one incident condition, θ_i = 40° and E_kin = 1.92 eV, is plotted. The black lines highlight a surface unit cell of α-Al₂O₃. Although the surface is mobile at the surface temperature of 300 K, the estimated location is accurate enough to analyse the origin of the bounce and the trajectories show a clear trend. In trajectories leading to small scattering angles the H-atom is mainly scattering from the topmost Al atoms, at medium angles the scattering predominately takes place at the O atoms and hollow sites between them, while at the largest scattering angles the scattering is again mainly occurring at the Al atoms. This scattering behavior can be explained by rather simple geometric considerations. In Fig. 8 a schematic representation of elastic scattering processes at different surface sites is shown for an incidence angle of θ_i = 40°. The top-layer Al atoms effectively shield the O layer preventing large and small scattering angles for the case of scattering at oxygen atoms, while collisions with the Al atoms can result in very low and large scattering angles. The rather flat oxygen layer thus only allows for scattering at angles close to the medium-sized incidence angle. These geometric conditions also apply to the other investigated incident conditions. Best agreement between simulation and experiment for E_kin = 1.92 eV is achieved for small and large scattering angels for which scattering from Al dominates. For intermediate scattering angles, where scattering from O dominates, the simulation results deviate more strongly from experiment. This might be a hint for a less accurate description of the scattering process at the oxygen sites. Since the comparison of the HDNNP trajectories with ab initio MD does show a similar good agreement for oxygen and aluminium sites, the most likely explanation thus is a less accurate description of the hydrogen–oxygen interaction by the RPBE functional.

Fig. 7 Impact positions of the H-atoms at the surface for the trajectories with incidence conditions θ_i = 40°, ϕ_i = 0° and an initial kinetic energy of 1.92 eV. For reference the equilibrated surface is shown. The black lines highlight the surface unit cell, darker grey circles are oxygen atoms, lighter grey circles are Al atoms. The brighter the color of the points, the higher is the density of impact sites.

Fig. 8 Model of the surface geometry resulting in different scattering angles for elastic scattering in case of an incidence polar angle of θ_i = 40°.

Having identified that trajectories initially scattered at oxygen atoms tend to show a too large energy loss, we took a closer look at representative trajectories to get a better understanding of its origin. Specifically, we analyzed the atomic kinetic energies and the changes of the positions of the surface atoms in selected trajectories by comparing them to trajectories with the same initial surface atom velocities but without the presence of a H-atom. In this way, we are able to see how all surface atoms are influenced by the collision with the H-atom. In this analysis we observed that if the H-atom hits an oxygen atom on the surface, the kinetic energy is transferred not only to one oxygen atom but is rapidly distributed over several neighboring surface atoms. Due to the involvement of multiple indirect collision partners during the collision with the O-atoms, the energy loss is increased compared to an idealized single-atom collision. Fig. S15 in the ESI† illustrates how the kinetic energy of the H-atom is efficiently transferred to multiple surface atoms. A similar mechanism has been observed for H-atom scattering from graphene.⁵² If, however, the H-atom is initially scattered at an Al-atom, the energy is mostly transferred to this Al with an overall much smaller energy loss in a more elastic process. This efficient energy transfer is also seen in ab initio MD simulations supporting that this effect is not a consequence of inaccuracies in the HDNNP.

One possible source of error contributing to deviations between experiment and simulation are surface imperfections. Although we took great care in experiment to prepare a well-defined surface, some defects are unavoidably present. From AFM measurements we know that the step density of the surface is about 2–4%. Furthermore, it is known that the α-Al₂O₃(0001) surface loses oxygen in an UHV environment. A significant O loss would lead to surface reconstruction and can be excluded, since we observed the (1 × 1) LEED pattern at all times. Still, oxygen vacancies will be present at the surface. Furthermore, the presence of hydrogen on the surface can not be excluded⁶⁰ and may affect the energy loss of the scattered hydrogen atoms. The deviation in the angular distributions is a clear indication that the surface is not perfectly flat. Steps and defects increase scattering normal to the surface, exactly what we observe (Fig. 6). However, our finding that the energy loss in experiment is smaller than in simulations can not easily be assigned to surface defects, as surface defects commonly lead to an increased energy loss, for example due to a higher probability of multi-bounce scattering. However, since a complex scattering mechanism involving many atoms is causing the larger energy loss at the oxygen sites in the simulations, it cannot be excluded that surface defects could change or even prohibit such a process.

The second point to be addressed is the general accuracy of the RPBE exchange correlation functional, which determines the interaction strength in between the surface atoms and between the surface atoms and the H atom. The best agreement between the experimental and theoretical kinetic energy distributions is obtained at low scattering angles with high kinetic energies. The largest difference in mean kinetic energy loss is present at large scattering angles and low kinetic incident energies, especially the scattering at E_kin,i = 0.99 eV and θ_s = 65°. Both cases represent extreme cases of H-atom-surface interaction. When the H-atom scatters at high kinetic energies and low scattering angles, the effective interaction time is minimal. However, when the scattering angle is large and the kinetic energy is low, the interaction time with the surface is significantly longer. As the kinetic energy distributions for high kinetic energies are in better agreement with the experiment compared to the distributions for low kinetic energies, apparently inaccuracies of the reference electronic structure method tend to accumulate with interaction time. Moreover, the bulk modulus of α-Al₂O₃ obtained in the DFT calculations is about 10% smaller than the experimental value (Table 1). This softer description of α-Al₂O₃ by the RBPE functional may lead to an overestimation of the kinetic energy loss.

5 Summary

We have studied the interaction of hydrogen atoms with the α-Al₂O₃(0001) surface using a combination of large-scale molecular dynamics simulations based on a high-dimensional neural network potential trained to DFT data employing the RPBE functional and highly accurate H-atom beam scattering experiments under UHV conditions. Two different initial kinetic energies and two different incidence polar angles have been investigated. Best agreement between experiment and theory is found for large initial kinetic energies and very low and high scattering angles, which we find to arise from scattering at top-layer aluminium atoms. Scattering at lower initial kinetic energies results in a larger loss of kinetic energy in the MD trajectories compared to experiment, and in general also scattering at oxygen sites seems to result in larger discrepancies between experiment and theory. As a careful validation of the multidimensional PES representation by the HDNNP showed a very close agreement with the underlying DFT reference calculations, we consider limitations in the accuracy of the employed RPBE functional as the most likely explanation of the observed deviations. Although, in view of the high kinetic energies, the differences between different exchange correlation functionals are small, more complex scattering mechanisms at the oxygen sites might increase the sensitivity of the simulations to subtle differences in the description of the atomic interactions. Moreover, future work should address the possible role of surface imperfections like steps, oxygen vacancies and adsorbed hydrogen atoms, which have not been included in the present work, on the energy loss and scattering mechanisms, which are more complex at the surfaces of binary compounds like α-Al₂O₃ compared to frequently studied elemental surfaces.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

O. B. is thankful to Alec M. Wodtke for continuous support in every aspect of the work. We thank Alexander Kandratsenka and Alec M. Wodtke for stimulating discussions. This project has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant numbers 389479699/GRK2455 BENCh and 217133147/SFB 1073, project A04 and C03. Open Access funding provided by the Max Planck Society.

Notes and references

1. Z. Ma and F. Zaera, Encyclopedia of Inorganic Chemistry – Heterogeneous Catalysis by Metals, John Wiley & Sons, Ltd, 2006, p. 1.
2. I. Fechete, Y. Wang and J. C. Védrine, Catal. Today, 2012, 189, 2.
3. J. G. Chen, R. M. Crooks, L. C. Seefeldt, K. L. Bren, R. M. Bullock, M. Y. Darensbourg, P. L. Holland, B. Hoffman, M. J. Janik, A. K. Jones, M. G. Kanatzidis, P. King, K. M. Lancaster, S. V. Lymar, P. Pfromm, W. F. Schneider and R. R. Schrock, Science, 2018, 360, eaar6611.
4. A. Tanimu, G. Tanimu, H. Alasiri and A. Aitani, Energy Fuels, 2022, 36, 5152.
5. L. Reichelt, S. Hippmann, V. N. Brichkin and M. Bertau, Z. Anorg. Allg. Chem., 2021, 647, 1583.
6. T. L. Jørgensen, H. Livbjerg and P. Glarborg, Chem. Eng. Sci., 2007, 62, 4496.
7. H. Chadwick and R. D. Beck, Annu. Rev. Phys. Chem., 2017, 68, 39–61.
8. G. B. Park, B. C. Krueger, D. Borodin, T. N. Kitsopoulos and A. M. Wodtke, Rep. Prog. Phys., 2019, 82, 096401.
9. O. Bünermann, A. Kandratsenka and A. M. Wodtke, J. Phys. Chem. A, 2021, 125, 3059.
10. G. Benedek and J. Toennies, Surf. Sci., 1994, 299, 587.
11. B. Gumhalter, Phys. Rep., 2001, 351, 1.
12. J. Werdecker, M. E. van Reijzen, B.-J. Chen and R. D. Beck, Phys. Rev. Res., 2018, 120, 053402.
13. Y. Alkoby, H. Chadwick, O. Godsi, H. Labiad, M. Bergin, J. T. Cantin, I. Litvin, T. Maniv and G. Alexandrowicz, Nat. Commun., 2020, 11, 3110.
14. A. Kleyn, A. Luntz and D. Auerbach, Phys. Rev. Res., 1981, 47, 1169.
15. T. Greenwood, H. AlSalem and S. P. K. Koehler, J. Phys. Chem. A, 2023, 127, 1124.
16. Z. Zhao, Y. Wang, X. Yang, J. Quan, B. C. Krueger, P. Stoicescu, R. Nieman, D. J. Auerbach, A. M. Wodtke, H. Guo and G. B. Park, Nat. Chem., 2023, 15, 1006.
17. L. Grajciar, C. J. Heard, A. A. Bondarenko, M. V. Polynski, J. Meeprasert, E. A. Pidko and P. Nachtigall, Chem. Soc. Rev., 2018, 47, 8307.
18. K. Reuter, C. Stampfl and M. Scheffler, Handbook of Materials Modeling, 2005, vol. 1, p. 149.
19. C. H. Christensen and J. K. Nørskov, J. Chem. Phys., 2008, 128, 182503.
20. A. Groß, Surf. Sci. Rep., 1998, 32, 291.
21. G. Wiesenekker, G. J. Kroes and E. J. Baerends, J. Chem. Phys., 1996, 104, 7344.
22. C. Crespos, M. A. Collins, E. Pijper and G. J. Kroes, Chem. Phys. Lett., 2003, 376, 566.
23. A. Groß, S. Wilke and M. Scheffler, Phys. Rev. Lett., 1995, 75, 2718.
24. C. M. Handley and P. L. A. Popelier, J. Phys. Chem. A, 2010, 114, 3371.
25. J. Behler, J. Chem. Phys., 2016, 145, 170901.
26. V. Botu, R. Batra, J. Chapman and R. Ramprasad, J. Phys. Chem. C, 2017, 121, 511.
27. V. L. Deringer, M. A. Caro and G. Csányi, Adv. Mater., 2019, 31, 1902765.
28. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K.-R. Müller, Chem. Rev., 2021, 121, 10142.
29. E. Kocer, T. W. Ko and J. Behler, Annu. Rev. Phys. Chem., 2022, 73, 163.
30. T. B. Blank, S. D. Brown, A. W. Calhoun and D. J. Doren, J. Chem. Phys., 1995, 103, 4129.
31. S. Lorenz, A. Groß and M. Scheffler, Chem. Phys. Lett., 2004, 395, 210.
32. J. Behler, S. Lorenz and K. Reuter, J. Chem. Phys., 2007, 127, 014705.
33. J. Behler, K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 115421.
34. I. Goikoetxea, J. Beltrán, J. Meyer, J. I. Juaristi, M. Alducin and K. Reuter, New J. Phys., 2012, 14, 013050.
35. J. Behler and M. Parrinello, Phys. Rev. Lett., 2007, 98, 146401.
36. J. Behler, Angew. Chem., Int. Ed., 2017, 56, 12828.
37. A. P. Bartók, M. C. Payne, R. Kondor and G. Csányi, Phys. Rev. Lett., 2010, 104, 136403.
38. A. V. Shapeev, Multiscale Model. Simul., 2016, 14, 1153.
39. J. S. Smith, O. Isayev and A. E. Roitberg, Chem. Sci., 2017, 8, 3192.
40. A. P. Thompson, L. P. Swiler, C. R. Trott, S. M. Foiles and G. J. Tucker, J. Comput. Phys., 2015, 285, 316.
41. O. T. Unke and M. Meuwly, J. Chem. Theory Comput., 2019, 15, 3678.
42. R. Drautz, Phys. Rev. B, 2019, 99, 014104.
43. S. Wille, H. Jiang, O. Bünermann, A. M. Wodtke, J. Behler and A. Kandratsenka, Phys. Chem. Chem. Phys., 2020, 22, 26113.
44. K. Shakouri, J. Behler, J. Meyer and G.-J. Kroes, J. Phys. Chem. Lett., 2017, 8, 2131.
45. P. Spiering, K. Shakouri, J. Behler, G.-J. Kroes and J. Meyer, J. Phys. Chem. Lett., 2019, 10, 2957.
46. N. Gerrits, K. Shakouri, J. Behler and G.-J. Kroes, J. Phys. Chem. Lett., 2019, 10, 1763.
47. C. L. Box, Y. Zhang, R. Yin, B. Jiang and R. J. Maurer, JACS Au, 2021, 1, 164.
48. O. Bünermann, H. Jiang, Y. Dorenkamp, A. Kandratsenka, S. M. Janke, D. J. Auerbach and A. M. Wodtke, Science, 2015, 350, 1346.
49. Y. Dorenkamp, H. Jiang, H. Kockert, N. Hertl, M. Kammler, S. M. Janke, A. Kandratsenka, A. M. Wodtke and O. Bunermann, J. Chem. Phys., 2018, 148, 034706.
50. B. C. Krueger, N. Bartels, C. Bartels, A. Kandratsenka, J. C. Tully, A. M. Wodtke and T. Schaefer, J. Phys. Chem. C, 2015, 119, 3268.
51. N. Hertl, A. Kandratsenka, O. Bünermann and A. M. Wodtke, J. Phys. Chem. A, 2021, 125, 5745.
52. H. Jiang, M. Kammler, F. Ding, Y. Dorenkamp, F. R. Manby, A. M. Wodtke, T. F. Miller, III, A. Kandratsenka and O. Buenermann, Science, 2019, 364, 379.
53. L. Lecroart, N. Hertl, Y. Dorenkamp, H. Jiang, T. N. Kitsopoulos, A. Kandratsenka, O. Buenermann and A. M. Wodtke, J. Chem. Phys., 2021, 155, 034702.
54. Y. Dorenkamp, C. Volkmann, V. Roddatis, S. Schneider, A. M. Wodtke and O. Bünermann, J. Phys. Chem. C, 2018, 122, 10096.
55. S.-H. Cho, N. D. Walther, S. T. Nguyen and J. T. Hupp, Chem. Commun., 2005, 5331.
56. T. Riittonen, E. Toukoniitty, D. K. Madnani, A.-R. Leino, K. Kordas, M. Szabo, A. Sapi, K. Arve, J. Wärnå and J.-P. Mikkola, Catalysts, 2012, 2, 68.
57. O. Bünermann, H. Jiang, Y. Dorenkamp, D. J. Auerbach and A. M. Wodtke, Rev. Sci. Instrum., 2018, 89, 094101.
58. L. Schnieder, K. Seekamprahn, F. Liedeker, H. Steuwe and K. Welge, Faraday Discuss., 1991, 91, 259.
59. C. C. Chang, J. Appl. Phys., 1968, 39, 5570.
60. J. Ahn and J. Rabalais, Surf. Sci., 1997, 388, 121.
61. G. Renaud, Surf. Sci. Rep., 1998, 32, 1.
62. D. P. Woodruff, Chem. Rev., 2013, 113, 3863.
63. C. Walters, K. McCarty, E. Soares and M. Van Hove, Surf. Sci., 2000, 464, L732.
64. J. Behler, Chem. Rev., 2021, 121, 10037.
65. J. Behler, J. Chem. Phys., 2011, 134, 074106.
66. J. Behler, Int. J. Quantum Chem., 2015, 115, 1032.
67. A. M. Tokita and J. Behler, J. Chem. Phys., 2023, 159, 121501.
68. J. Behler, J. Condens. Matter Phys., 2014, 26, 183001.
69. V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter and M. Scheffler, Comput. Phys. Commun., 2009, 180, 2175.
70. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
71. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
72. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413.
73. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158.
74. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
75. M. Herbold and J. Behler, J. Chem. Phys., 2022, 156, 114106.
76. N. Artrith and J. Behler, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 45439.
77. J. S. Smith, B. Nebgen, N. Lubbers, O. Isayev and A. E. Roitberg, J. Chem. Phys., 2018, 148, 241733.
78. E. V. Podryabinkin and A. V. Shapeev, Comput. Mater. Sci., 2017, 140, 171–180.
79. Y. Eldar, M. Lindenbaum, M. Porat and Y. Zeevi, IEEE Trans. Image Process., 1997, 6, 1305.
80. T. B. Blank and S. D. Brown, J. Chemom., 1994, 8, 391–407.
81. S. Plimpton, J. Comput. Phys., 1995, 117, 1.
82. A. Singraber, T. Morawietz, J. Behler and C. Dellago, J. Chem. Theory Comput., 2019, 15, 3075.
83. W. C. Swope, H. C. Andersen, P. H. Berens and K. R. Wilson, J. Chem. Phys., 1982, 76, 637.
84. D. J. Evans and B. L. Holian, J. Chem. Phys., 1985, 83, 4069.
85. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1.
86. A. Togo, Phys. Soc. Jpn., 2023, 92, 012001.
87. X.-L. Wang, C. R. Hubbard, K. B. Alexander, P. F. Becher, J. A. Fernandez-Baca and S. Spooner, J. Am. Ceram. Soc., 1994, 77, 1569.
88. D. L. Anderson and O. L. Anderson, J. Geophys. Res., 1970, 75, 3494.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cp04729f

---