Best-of-both-worlds computational approaches to difficult-to-model dissociation reactions on metal surfaces

Abstract

The accurate modeling of dissociative chemisorption of molecules on metal surfaces presents an exciting scientific challenge to theorists, and is practically relevant to modeling heterogeneously catalyzed reactive processes in computational catalysis. The first important scientific challenge in the field is that accurate barriers for dissociative chemisorption are not yet available from first principles methods. For systems that are not prone to charge transfer (for which the difference between the work function of the surface and the electron affinity of the molecule is larger than 7 eV) this problem can be circumvented: chemically accurate barrier heights can be extracted with a semi-empirical version of density functional theory (DFT). However, a second important challenge is posed by systems that are prone to (full or partial) electron transfer from the surface to the molecule. For these systems the Born–Oppenheimer approximation breaks down, and currently no method of established accuracy exists for modeling the resulting effect of non-adiabatic energy dissipation on the dissociative chemisorption reaction. Because two problems exist for this class of reactions, a semi-empirical approach to computing barrier heights, which would demand that computed and experimental dissociative chemisorption probabilities match, is unlikely to work. This Perspective presents a vision on how these two problems may be solved. We suggest an approach in which parameterized density functionals are used as in the previous semi-empirical approach to DFT, but in which the parameters are based on calculations with first principles electronic structure methods. We also suggest that the diffusion Monte-Carlo (DMC) and the random phase approximation (RPA) probably are the best two first principles electronic structure methods to pursue in the framework of the approach that we call first-principles based DFT (FPB-DFT) – providing DMC and the RPA with a steppingstone towards benchmarking and future applications in computational catalysis. Probably the FPB density functional is best based on screened hybrid exchange in combination with non-local van der Waals correlation. We also propose a new electronic friction method called scattering potential friction (SPF) that could combine the advantages and avoid the disadvantages of the two main existing electronic friction approaches for describing non-adiabatic effects: by extracting an electronic scattering potential from a DFT calculation for the full molecule–metal surface system, it might be possible to compute friction coefficients from scattering phase shifts in a computationally convenient and robust fashion. Combining the FPB-DFT and SPF methods may eventually result in barrier heights of chemical accuracy for the difficult-to-model class of systems that are prone to charge transfer. This should also enable the construction of a representative database of barrier heights for dissociative chemisorption on metal surfaces. Such a database would allow testing new density functionals, or, more generally, new electronic structure approaches on a class of reactions that is of huge importance to the chemical industry. Additionally, the difficult-to-model sub-class of systems we focus on is essential to sustainable chemistry and important for a sustainable future. Adding the database envisaged to large databases already existing but mostly addressing gas phase chemistry will enable testing density functionals that have a claim to universality, i.e., to be good for all chemical systems of importance. We also make a suggestion for how to develop such a generally applicable functional, which should have the correct asymptotic dependence of the exchange contribution to the energy in both the gas phase and the metal. Finally we suggest some improvements in the representation of potential energy surfaces and in dynamics methods that would help with the validation of the proposed methods.

---

Geert-Jan Kroes

Geert-Jan Kroes obtained his PhD in Chemistry in 1990 at the University of Amsterdam (The Netherlands), working under the supervision of Prof. R. P. H. Rettschnick. He worked as a post-doc with David Clary from 1990 to 1992 (Cambridge, UK), and with Marc van Hemert and Ewine van Dishoeck from 1992 to 1993 (Leiden, The Netherlands). He worked as a KNAW-fellow at the Free University of Amsterdam and at Leiden University from 1993–1998. At Leiden he became Assistant Professor in 1998, and Full Professor in 2003. His research is focused on reactive scattering of molecules from metal surfaces.

Jörg Meyer

Jörg Meyer's main research interest is to understand energy dissipation at the atomic scale, particularly in the context of chemical reactions. Surfaces provide a variety of energy dissipation channels, which are constituted by the nuclear and electronic degrees of freedom of the surface atoms. He uses a variety of computational modelling techniques, from electronic structure to force field levels, and investigates important implications, for example, in heterogeneous catalysis and astrochemistry.

---

1. Introduction

The production of most chemicals involves heterogeneous catalysis at some stage.¹ Heterogeneously catalyzed processes usually involve a sequence of elementary reactions on surfaces, and often these surfaces are of metal particles. An important class of elementary reactions is formed by dissociative chemisorption reactions, in which a bond in the molecule breaks and both resulting fragments bond to the surface. Such reactions are interesting for both practical and scientific reasons.

Minimum barrier heights (E_b) (or transition states, TSs) for dissociative chemisorption (DC) are not just important for this specific reaction alone. Heterogeneously catalysed reactions consist of a sequence of elementary reaction steps. Among these sequences barriers to DC reactions often control the rate (we should actually say: the “turn-over-frequency”) of an entire heterogeneously catalysed process.^2,3 Important examples are steam reforming for hydrogen production⁴ and production of ammonia to make fertilizer.^5,6 The accuracy with which these E_b (TSs) can be calculated has an important impact on computational heterogeneous catalysis.^7,8 With the present accuracy one can predict trends in heterogeneous catalysis using transition metals (TMs), and, with reasonable accuracy, which catalysts should be good for producing specific chemicals.^8,9 But theory is not yet very accurate for turn-over-frequencies of entire heterogeneously catalysed processes, as illustrated by errors in the computed rate of ammonia production still being as large as 1–2 orders of magnitude.^8,10 This poses a major problem to the accurate calculation of these rates, and the catalysis literature has emphasized the need for higher accuracy than afforded now by current density functionals (DFs) for achieving predictive power.¹¹ It is hard to overemphasize the importance of this. Catalysis is “the single most important and pervasive interdisciplinary technology in the chemical industry”, it has been a “key enabling technology” for most of the societal challenges in Horizon 2020, and with the production of artificial fertilizers HC makes it possible to feed a world population of 7 billion humans.¹² Improvements in catalysts can reduce the energy demand of the chemical industry by at least 20–40%.¹²

The study of how molecules react at or, more generally, scatter from metal surfaces is also to a large extent curiosity driven. The study of these processes is the subject of the field of molecule–surface reaction dynamics.^13–19 Instead of reacting a molecule may also scatter back to the gas phase, possibly changing its vibrational and/or rotational state. These outcomes all depend on e.g. the molecule's incident translational energy, its initial rovibrational state, the surface temperature (T_s), and, crucially, on how the molecule interacts with the surface.¹³ Whether or not the molecule will react will of course depend on the barrier height to reaction but, just like for gas phase reactions,^20–25 the accurate prediction of barriers for DC reactions on metals comes with major challenges.^13,26–28 The latter reactions have been addressed mostly with density functional theory (DFT), which is perhaps not surprising: these reactions take place in complex systems containing many electrons to be modelled, and with more than 30 000 papers published annually DFT is probably the most important electronic structure method for complex systems.²⁹ What is surprising, though, is that far more databases with benchmark data are available for testing new electronic structure methods on gas phase reaction barriers^23–25,30 than on dissociative-chemisorption-on-metals barriers.^31,32 Given the large practical importance of the latter category of reactions, one might expect the opposite. As will become clear from the below, this has to do with the fact that accurate quantum chemical calculations are much more challenging to carry out for molecules interacting with metal surfaces than for molecules interacting with one another in the gas phase. Before we dive into explaining this, let us first explain what the present state-of-the-art is in the former area of applications.

To determine barrier heights (E_b) for DC on metals, in the present state-of-the-art a semi-empirical electronic structure approach, called specific reaction parameter DFT (hence: SRP-DFT),^13,33–35 is taken. Here one takes a judiciously chosen^13,36 density functional (DF) with a single adjustable parameter. This parameter is then fitted so that calculations based on the adjusted DF and using a suitable dynamics model and method reproduce DC probabilities that are measured as a function of translational incidence energy in supersonic molecular beam experiments. In Surface Science DC probabilities are also called by the shorter name “sticking probabilities” and denoted by the symbol S₀, as we will often do from now on. Even though only one parameter, which correlates with E_b, in the DF is tweaked in SRP-DFT, the approach yields DFs that are generally capable of reproducing S₀ curves over large ranges of energies. The reason is that standard DFT, though not capable of predicting E_b for a DC-on-metal system accurately, is capable of yielding quite a good description of the variation of the barrier height with the system's configuration (e.g., the molecule's orientation and impact site on the surface, for a diatomic molecule).^37,38 According to the so-called hole model³⁹ this should be a sufficient condition for theory to be able to compute accurate S₀, although in practice it is also necessary to choose a dynamical model and method capable of accurately describing dynamical effects.^13,40 Using the SRP-DFT approach, it has been possible to extract E_b for 14 systems in which a molecule dissociatively chemisorbs on a metal surface, which has led to a database that can be used to test existing and new DFs on their accuracy for DC for these systems.³²

So far, in the SRP-DFT approach the exchange part of the exchange–correlation functional of DFT has been based on the generalized gradient approximation (GGA).²⁰ Unfortunately, this approach breaks down for systems that are characterized by the charge transfer energy E(CT) being smaller than 7 eV.⁴¹ Here, E(CT) is defined as the work function of the metal surface minus the electron affinity of the incident molecule. In systems with E(CT) < 7 eV the metal easily lets go of an electron (low work function) and/or the molecule is eager to accept it (a positive electron affinity implies a high affinity for electrons). As a result such systems are prone to (partial) electron transfer from the surface to the molecule. GGA functionals perform poorly for systems affected by charge transfer.⁴² Such breakdown is illustrated for an individual system in Fig. 1, for HCl + Au(111), which has E(CT) = 5.8 eV. Even using about the most repulsive GGA DF still compatible with a good description of the metal (RPBE⁴⁸), the measured DC probability is overestimated, suggesting that E_b is underestimated with the RPBE DF. This holds true even though RPBE typically overestimates E_b for systems with E(CT) ≥ 7 eV.³² Recent work⁴⁹ has suggested that the difference might be due to the calculations using a higher “coverage” than the experiments, which were done in the zero-coverage limit.⁴³ We believe that this is not correct: in the recent experiments, “coverage” refers to pre-coverage by Cl.⁴⁹ In the calculations,⁴³ “coverage” refers to the density of the incident HCl beam. AIMD calculations found little difference between sticking probabilities employing a (3 × 3) surface unit cell and a (2 × 2) surface unit cell (HCl “coverages” of 1/9 and 1/4, respectively).⁴⁵ This result implies that the (2 × 2) dynamics calculations should already be useable in comparisons with experiments in the zero-coverage limit.

Fig. 1 Sticking probabilities computed⁴³ for HCl + Au(111) with the MS-RPBEl meta-GGA DF,⁴³ the RPBE GGA DF,⁴⁴ and the SRP32-vdW1 DF featuring semi-local exchange and non-local correlation⁴⁵ are compared to experimental results, for normal incidence. The open green squares are the experimental sticking probabilities first published,⁴⁶ while the upper base and lower base triangles represent upper and lower bounds to the experimental S₀ obtained from an improved analysis of the experiments.⁴³ Panel (a) uses a linear and panel (b) a logarithmic scale for S₀. In panel (b) results of QCT⁴⁴ (orange diamonds connected by solid orange line) as well as of QD⁴⁷ (orange dashed lines) calculations using the RPBE DF are presented. Figure taken from ref. 43 (https://doi.org/10.1021/acs.jpcc.0c03756). Further permission requests to be directed to the ACS.

As illustrated by Fig. 2, similar observations as for HCl + Au(111) have been made for other systems with, as a rule of thumb, E(CT) < 7 eV,⁴¹ such as O₂ + Al(111)^51,52 and Cu(111),⁵³ H₂O + Ni(111),⁵⁴ and NH₃ + Ru(0001).^55,56 By the phrase “rule of thumb” we mean to say that E(CT) 7 eV can be taken as an approximate number for the value of E(CT) that determines whether a GGA likely breaks down according to the empirical evidence presented in Fig. 2 and ref. 41. In line with this trend,⁴¹ RPBE-DFT calculations underestimated the barrier height computed for H₂ + Mg(0001) with diffusion Monte-Carlo,⁵⁷ which is an accurate stochastic first-principles-based electronic structure method.⁵⁸ As also illustrated by Fig. 2 for systems with E(CT) > 7 eV it has been possible to model the DC on metals with chemical accuracy using GGA density functionals.

Fig. 2 Correlation of the ability of DFs based on GGA-exchange with E(CT), i.e., the difference between the work function of the metal surface and the electron affinity of the incoming molecule⁵⁰ (shown with vertical lines). The blue and green lines represent systems for which it was possible to derive a SRP DF. The red lines represent systems for which the use of the repulsive RPBE GGA DF leads to overestimating computed S₀, while the orange lines indicate systems for which computed results strongly suggest that this is the case. See Fig. 70 of ref. 13 for how the work functions of the surfaces and the electron affinities of the molecules were obtained. Figure taken from ref. 41 (https://pubs.acs.org/doi/10.1021/acs.jpclett.0c02452). Further permission requests to be directed to the ACS.

As noted in ref. 41 and in the caption of Fig. 2 the evidence for GGA-breakdown for E(CT) < 7 eV is solid for some systems (indicated with red in Fig. 2) while still subject to some uncertainty for other systems (indicated with orange in Fig. 2). For instance, in quantum dynamics calculations on D₂O + Ni(111) making approximations to parallel translational motion of the molecule and the surface vibrations and based on a RPBE⁴⁸ GGA potential energy surface, the RPBE results overestimated experimental sticking probabilities for the 2ν₃ state, while they underestimated the measured values for the 1ν₃ state at low average incidence energies (see Fig. 6 of ref. 54). It is not clear which comparison should be given precedence: the experimental probabilities are larger for the 2ν₃ state (at values up to about 10⁻² instead of values up to about 10⁻³ for the 1ν₃ state), suggesting that the comparison should be the most important for the higher excited vibrational state. On the other hand errors in the shape of the PES, which could also occur for a GGA functional, might make it harder to accurately describe the reactivity of higher excited vibrational states with theory. In conclusion, the “border” of 7 eV should not be interpreted strictly, but rather as a rule of thumb. More accurate dynamical calculations (i.e., using less dynamical approximations) may be required than was possible in 2016 (ref. 54) to determine whether RPBE fails for D₂O + Ni(111), for which E(CT) ≈ 5.4 eV. Note that an accurate comparison cannot be based on experimental “laser-off” results and theoretical results for the vibrational ground state, as an appropriate inclusion of excited vibrational states would lead to a theoretical laser-off reactivity that should be much higher than the theoretical result for the ground vibrational state only.^54,59

A similar problem as discussed above for systems with E(CT) < 7 eV is observed when modeling experiments on vibrationally inelastic scattering, especially if the incident molecule is initially in a highly excited vibrational state. Vibrationally inelastic scattering, in which the molecule's vibrational state changes in an otherwise non-reactive molecule–surface collision, might be thought to lie in the realm of physics rather than chemistry. The process is however intimately related to DC, as changes in the bond length of the incident molecule are implicated in both processes and thus probe the barrier region of the potential energy surface (PES) of the system. This interrelatedness is also illustrated in Fig. 3, which compares computational with experimental results for vibrationally inelastic scattering of NO from Au(111), for which E(CT) = 5.3 eV.⁶³ In Fig. 3, the calculations labeled BOMD, MDEF(LDFA), MDEF(ODF) used a PES that allowed for DC, i.e., bond stretching can result in reaction.⁶⁴ In contrast, the calculations labeled MDEF and IESH used a PES that did not allow reaction.⁶⁵ Clearly the use of these two different PESs leads to different results, but Fig. 3 also illustrates another problem for systems with E(CT) < 7 eV. As mentioned earlier, in these systems the metal tends to easily let go of an electron (low work function) and the molecule is eager to accept it (again, a positive electron affinity implies a high affinity for electrons). This implies that the system is prone to non-adiabatic effects related to electron–hole pair (ehp) excitation, as also illustrated in Fig. 4.¹⁵ In Fig. 3 the acronyms used are telltale of how these non-adiabatic effects are described. “MDEF” stands for molecular dynamics calculations describing ehp excitation in a weak coupling approximation using “electronic friction” (EF).⁶⁶ “IESH” calculations used the independent electron surface hopping (IESH) method, which allows for strong coupling between two diabatic electronic states of the molecule–metal system (i.e., neutral-metal and anion–metal) with ehp excited states superimposed.^67,68 In spite of the improvements made to the PES,⁶⁴Fig. 3 therefore also illustrates how the present state-of-the-art theory dealing with non-adiabatic effects falls short in describing the chemical physics problem of vibrationally inelastic scattering of NO from Au(111).

Fig. 3 Experimental^60–62 relative probabilities for scattering of NO from Au(111) in its initial vibrational state v_i to the final vibrational state v_f are compared with calculations⁶³ using different dynamical models and potential energy surfaces. Results are given for (a) v_i = 2, (b) v_i = 3, (c) v_i = 11, and (d) v_i = 16. The reference in the legend corresponds to the reference with that number in ref. 63. The relevant details are discussed in the text, for other details see ref. 63. Figure taken from ref. 63 (https://doi.org/10.1021/jacsau.0c00066). Further permission requests to be directed to the ACS.

Fig. 4 Correlation of E(CT) with non-adiabatic effects. (a) Non-adiabatic effects are governed by E(CT), i.e., the difference of the work function of the surface, here indicated with the symbol Φ, and the electron affinity of the molecule, EA, as indicated. The energy needed to transfer an electron from the surface to the molecule can be diminished by the image–charge interaction, which stabilizes the indicated anion affinity level of the approaching molecule. If needed the remaining energy required for full electron transfer can come from e.g. the incident translational and vibrational energy of the impinging molecule. (b) E(CT) as an indicator of how likely non-adiabatic effects are. Systems to the left (small E(CT)) are more prone to ehp excitation than systems to the right. For details see ref. 15. Used with permission of Annual Reviews of Physical Chemistry, from [The Dynamics of Molecular Interactions and Chemical Reactions at Metal Surfaces: Testing the Foundations of Theory, K. Golibrzuch, N. Bartels, D. J. Auerbach, and A. M. Wodtke, vol. 66, pp. 399–425, 2015]; permission conveyed through Copyright Clearance Center, Inc.

The problem with describing electronically non-adiabatic effects can be illustrated for DC of N₂ on Ru(0001),⁶⁹ for which E(CT) at 7.4 eV actually somewhat exceeds 7 eV, as illustrated in Fig. 5. As in Fig. 2 “LDFA and “ODF” label two different EF methods that can be used to describe ehp excitation; below we will outline their differences and describe them in some detail. The point made by Fig. 5 is that the computed S₀, which are shown on a log scale, depend dramatically on which of the two EF methods (LDFA or ODF) is used. Specifically, the S₀ curves computed with molecular dynamics with electronic friction (MDEF) using the LDFA and the ODF are removed from one another in incidence energy by more than 2 kcal mol⁻¹ (1 kcal mol⁻¹ ≈ 0.043 eV). This is problematic because at present it is unclear if any on the two EF methods is capable of accurately describing the effect of ehp excitation on DC on metals.

Fig. 5 S ₀ computed for N₂ + Ru(0001) with MDEF and the RPBE DF using the Born–Oppenheimer static surface (BOSS) model,⁷⁰ the Born–Oppenheimer moving surface (BOMS) model,⁷⁰ and the non-Born–Oppenheimer moving surface (NBOMS) model with the LDFA and the ODF approximations⁶⁹ are shown as a function of E_i. Experimental results are also shown. Figure taken from ref. 69 (https://pubs.acs.org/doi/10.1021/acs.jpclett.9b00523). Further permission requests to be directed to the ACS.

The above observations lead to the following clear conclusion for systems with E(CT) < 7 eV: given the uncertainties surrounding the accuracy of the present methods developed to describe non-adiabatic effects, an approach to deriving the PES with a semi-empirical DF based on reproducing measured sticking probabilities alone is clearly doomed to fail (see below). In the present state-of-the-art the choice of the semi-empirical DF is likely to depend on what method is used to describe non-adiabatic effects, which are strong for E(CT) < 7 eV (Fig. 4). At the same time, as will be discussed below first principles electronic structure methods with demonstrated chemical accuracy for dissociative-chemisorption-on-metal systems in general are not yet available, and the workhorse method (DFT-GGA) fails for systems with E(CT) < 7 eV (Fig. 2). Therefore, a risk exists for such systems that if agreement with experiment is achieved the agreement might be achieved “for the wrong reason”, i.e., due to cancellation of errors resulting from errors in both the PES and in the method used to describe non-adiabatic effects. In any case for such systems a semi-empirical approach is bound to go south if GGA-DFT is used, in which case “delocalization” errors are expected in the DFT due to (partial) charge transfer,⁴² which additionally causes strong non-adiabatic effects¹⁵ that we do not yet know how to model accurately.⁶⁹ In other words, when E(CT) < 7 eV it never rains but it pours.

At this point it is good to stop for a moment and to also emphasize the practical importance of systems for which E(CT) < 7 eV. It will be clear that these systems will often involve molecules with a high affinity for electrons, including oxygen containing molecules like O₂, methanol, and CO₂. It is hard to overemphasize the practical importance of these systems. The oxygen reduction reaction is crucial to clean energy conversion in fuel cells, and the chemistry of O₂ is of tremendous overall importance.⁷¹ Methanol and related compounds such as dimethyl ether and oxy-methylene ethers, which can be produced from CO₂,⁷² are considered to be important clean energy fuels in the future.¹² In short, many of the systems that pose the challenges discussed above are also crucial to sustainable chemistry, and it is therefore also practically important to address these challenges.

As discussed below, the greatest challenge is to come up with an accurate electronic structure approach that is ultimately based on first principles, and which can be used for predictive purposes. As emphasized elsewhere³⁷ such a method is also preferred in general for other reasons: standing in the way of semi-empirical approaches, experiments might yield conflicting results, they may not be well documented, or they may simply be absent for the system of interest.¹³ More accurate PESs will also be needed if the field studying collisions of molecules with metal surfaces is ever to match the field studying gas phase molecular collisions for level of detail and accuracy.⁷³ And once it will be possible to accurately compute PESs it will also be possible to test different methods for dealing with non-adiabatic effects more rigorously than now feasible for systems with E(CT) smaller than or approximately equal to 7 eV. Developing improved methods for dealing with such non-adiabatic effects also represents a very important challenge, which can already be addressed now. Finally, a number of challenges is related to how dynamics calculations need to be performed to validate new theories through comparison of computed observables with experiments.

All of the challenges described above will be addressed below. For this we will first address the state-of-the-art in electronic structure theory (Section 2.1). Next, we will consider theories dealing with electronically non-adiabatic effects in scattering of molecules from metal surfaces (Section 2.2). After that we will briefly consider methods needed to compute sticking probabilities with dynamics, including methods for PES construction (Section 2.3). Next, we will sketch a vision of how to proceed with improvements. Improvements in electronic structure theory and in methods for dealing with non-adiabatic effects will be addressed in Sections 3.1 and 3.2. Improvements in the representation of PESs and in dynamics calculations will be briefly discussed in Sections 3.3. Section 3.4 deals with the validation of the new methods through comparison with experiments. We end with summary conclusions (Section 4).

2. The state-of-the-art in modeling dissociative chemisorption on metals

2.1 State-of-the-art electronic structure theory

Present theoretical research on DC of molecules on metal surfaces usually bases itself on electronic structure calculations performed with DFT, using DFs at the GGA rung of Jacob's ladder,²⁰ and in some cases employing DFs at the meta-GGA rung,²⁰ both of which belong to the group of semi-local functionals. GGA and meta-GGA DFs are called semi-local because they evaluate the exchange–correlation energy density at a point in space as a function of that electronic density and its derivatives.⁷⁴ In some cases GGA correlation or meta-GGA correlation is replaced with a non-local correlation functional. For this typically the vdW-DF1 (ref. 75) or vdW-DF2 (ref. 76) correlation functionals are combined with GGA exchange, or the rVV10 (ref. 77) correlation functional is combined with meta-GGA exchange. The problem usually encountered with such functionals is illustrated for strictly semi-local functionals (i.e., also using semi-local correlation) in Fig. 6. This figure shows results of dynamical calculations of the probability of sticking in the benchmark system H₂ + Cu(111)⁷⁹ with comparison to experiment.⁷⁸ According to the hole model,³⁹ which states that for a given incidence energy the sticking probability should equal the proportion of geometries (impact site and molecular orientation) for which the molecule's energy exceeds the geometry's barrier height, agreement between theory and experiment of S₀ would signal the accurate calculation of barriers. As will be clear from Fig. 6 the accurate computation of barriers is accomplished with none of the GGA (PBE⁸⁰ and RPBE⁴⁸) and meta-GGA (TPSS,⁸² revTPSS,⁸³ and SCAN⁸¹) standard functionals tested. We also use Fig. 6 to reiterate a point already alluded to implicitly above: barrier heights are not observables. Their validation requires dynamics calculations of measured S₀ as illustrated with Fig. 6.

Fig. 6 Values of S₀ measured⁷⁸ in molecular beam experiments for D₂ + Cu(111) are compared to computed⁷⁹ values using PESs based on DFT calculations using the GGA PBE⁸⁰ and RPBE⁴⁸ DFs, and the meta-GGA SCAN,⁸¹ TPSS,⁸² and rev-TPSS⁸³ DFs. Figure taken from ref. 79 (https://doi.org/10.1021/acs.jpca.9b02914). Further permission requests to be directed to the ACS.

As alluded to in the Introduction the type of standard functionals discussed above can be combined in a semi-empirical way to model a particular (“easy-to-model”, E(CT > 7 eV) class of systems with chemical accuracy. A generic form of a parameterizable exchange–correlation functional that can be used for this is

E_XC = xE_X1 + (1 − x)E_X2 + E_C (1)

Here, E_Xi denotes an exchange functional, E_C a correlation functional, and E^DF_X and E^DF_C are exchange and correlation parts of the exchange–correlation functional DF. In writing eqn (1) we have assumed that the exchange–correlation functional can be separated in an exchange and a correlation part. Typical choices^13,36 that can be made for the exchange functionals in the semi-empirical SRP approach to DFT are E_X1 = E^RPBE_X, E_X2 = E^PBE_X or E^PBEsol_X,⁸⁴ and E_C = E^PBE_C, E^vdW-DF1_C,⁷⁵ or E^vdW-DF2_C,⁷⁶ but other choices are possible.^32,79,85Fig. 7 presents a typical example (D₂ + Pt(111))⁸⁶ showing that for easy-to-model systems agreement to within chemical accuracy (1 kcal mol⁻¹ ≈ 4.2 kJ mol⁻¹) can be obtained with experiment for the sticking probability computed using a PES calculated with SRP-DFT.

Fig. 7 Reaction probabilities computed⁸⁶ with SRP-DFT for D₂ + Pt(111) are compared with experimental values.⁸⁷ The green symbols indicate at which interpolated collision energies measured sticking probabilities would be obtained for values equal to those computed with theory. The numbers indicate the distance (in kJ mol⁻¹) of the corresponding points on the interpolated experimental curve to the points computed with dynamics along the collision energy axis. Reprinted from E. N. Ghassemi, M. Wijzenbroek, M. F. Somers, and G. J. Kroes, Chemically accurate simulation of dissociative chemisorption of D₂ on Pt(111), Chem. Phys. Lett., 2017, 683, 329–335, licensed under CC-BY 4.0.

With the SRP-DFT procedure a database (called SBH17) has been constructed presenting accurate barrier heights for DC in 17 molecule–metal surface systems, and this database has been used to benchmark the performance of 14 density functionals.³² We show results for a selection of these functionals (2 GGAs, 2 meta-GGAs, and 2 combinations of GGA exchange with non-local van der Waals correlation) in Table 1. Notable conclusions from Table 1 are that the best performing DF was the general purpose PBE⁸⁰ GGA density functional, and that the MS2 “made-simple” DF⁸⁸ was the best performing meta-GGA DF, perhaps because it was explicitly constructed to do well for metals as well as molecules. Importantly, the RPBE⁴⁸ and BEEF-vdW2 (ref. 89) functionals do not perform so well on barriers of easy-to-model systems. The performance of these functionals on barriers has been overrated due to systematic errors made with computing barrier heights in an earlier study that resulted in the SBH10 database.³¹ In this SBH10 study activation energies computed with DFT were compared with reference values of classical barrier heights for 9 of the 10 systems selected.³¹ There are two systematic differences between how the reference values of classical barrier heights and the values of the activation energies have been computed:³¹ (i) the activation energies contain zero-point energy corrections for the reactants and the transition state (barrier geometry), whereas the classical barrier heights do not. (ii) In the calculations of the activation energies the surface atoms were allowed to relax in the presence of the dissociating molecule in the transition state.³¹ This was not allowed in the SRP-DFT calculations of the reference values of classical barrier heights.³² The reason is that, so far, SRP-DFT has been built on the comparison with sticking experiments using supersonic molecular beams, and in the underlying dynamics of DC the surface atoms do not have time to relax to the full minimum barrier geometry (in which the surface atoms would be relaxed).¹³ This is because the surface atoms usually do not have time to respond to the (fast) incoming molecule in a supersonic molecular beam experiment.

Table 1 Performance for a selection of 6 density functionals (out of a batch of 14 functionals tested) on the SBH17 database of dissociative chemisorption barriers on metals.³² The MAE and MSE are in kcal mol⁻¹. The r_MAE and the r_MSE are the ranks describing how well the DFs performed out of the 14 DFs originally tested for these error measures, with a rank of 1 indicating best performance. For details see ref. 32

Functional	Type of DF	MAE	r MAE	MSE	r MSE
PBE	GGA	2.38	1	−1.34	5
RPBE	GGA	5.26	13	5.26	13
PBEα57-vdW-DF2	GGA + vdW	2.86	4	−0.92	3
BEEF-vdW-DF2	GGA + vdW	4.40	10	4.40	10
MS2	Meta-GGA	2.70	3	−1.71	6
SCAN	Meta-GGA	3.23	7	−2.42	8

---

Fig. 8 illustrates both a reason for why the SRP procedure based on semi-local exchange DFs works so well for easy-to-model systems (E(CT) > 7 eV), and why it will usually fail for systems that are difficult-to-model (for which, as a rule of thumb, E(CT) > 7 eV). Fig. 8 presents results⁵⁸ for an easy-to-model system (H₂ + Cu(111), E(CT) = 8.1 eV) and a system that turned out to be difficult-to-model (H₂ + Al(110), E(CT) = 7.2 eV). For H₂ + Cu(111) barrier heights computed with GGAs, meta-GGAs, and functionals containing GGA exchange and non-local correlation scatter around the semi-empirical SRP-DFT value. It thus makes sense that appropriately mixing such functionals according to eqn (1) allows an experiment on sticking of H₂ on Cu(111) to be reproduced. For H₂ + Al(110), on the other hand, none of these functionals, including the more repulsive ones, overestimates the semi-empirical value of the barrier height.⁵⁸ It is for the same reason that DC of HCl on Au(111) has, so far, not been well described (see Fig. 1): the functionals used up to now all contained semi-local exchange and therefore underestimated barrier heights resulting in overestimated computed S₀. Clearly, for difficult-to-model systems an approach that is more generally applicable than combining semi-local exchange with semi-local or non-local correlation from DFT is needed.

Fig. 8 Lowest reaction barrier heights for H₂ + Al(110) vs. those for H₂ + Cu(111), both taken relative to the best semi-empirical values, as obtained with DFs including semi-local exchange (i.e., not including exact exchange, indicated by blue symbols), with screened hybrid DFs (yellow and orange symbols), with diffusion Monte-Carlo (red symbol) and with the RPA (maroon symbol). The concentric gray shaded ellipses indicate areas that would be circles in an equal-distanced representation, with radii of 1 (darkest grey), 2, 3, and 4 (lightest gray) kcal mol⁻¹ respectively. For details see ref. 58. Figure reprinted from B. Oudot, K. Doblhoff-Dier, Journal of Chemical Physics, vol. 161, Article ID 054708, 2024; licensed under a Creative Commons Attribution (CC-BY) license.

For gas phase reactions semi-local functionals systematically underestimate barrier heights.^20–25 This has variously been attributed to self-interaction errors and delocalization errors,^20–22,42 where the latter may be considered an overarching term encompassing self-interaction error.⁴² To improve the description of gas phase barriers a successful strategy has been to use hybrid functionals, which are on the fourth rung of the DFT ladder²⁰ and admix exact exchange (EXX, also called Hartree–Fock exchange) with semi-local exchange.^23–25 In this endeavor, it has often paid off to use a somewhat higher fraction of exact exchange than used in hybrid functionals for general purposes.^23,42,90,91 When dealing with metals it is important to screen the EXX at long distances between the electrons, to avoid a collapse of the density of states of the electrons at the Fermi-level.⁹² For DC of O₂ on Al(111), which is an infamously difficult-to-model system, it has recently been shown⁴¹ that sticking probabilities can be obtained with a much higher accuracy than possible with semi-local DFs if a screened hybrid DF with an increased fraction of EXX (HSE03-1/3X)^41,93,94 is used, as shown in Fig. 9. While it is yet to be demonstrated that screened hybrids also perform systematically better for easy-to-model systems, it is encouraging to see that the very similar screened hybrid HSE06-1/3X^58,97 outperformed all GGAs and meta-GGAs tested on the two H₂ + Al(110) and Cu(111) systems (see also Fig. 8). The observation to the contrary in tests on the SBH10 database³¹ can most likely⁵⁸ be attributed to the incorrect comparison of computed activation energies with reference values of classical barrier heights, as detailed before. While calculations using hybrid functionals are expensive even with screened hybrids, it is also encouraging to note that at least for O₂ + Al(111) the use of a non-self-consistent approach using GGA densities (HSE03-1/3X) was quite successful.⁹⁸ We note that electronically adiabatic dynamics calculations were used to compare dynamics results based on the HSE03-1/3X PES with experiment, thus neglecting ehp excitation. The reasons that we expect this approximation to be accurate for O₂ + Al(111) is that good agreement with experiment was also obtained with adiabatic dynamics calculations using an accurate PES based on embedded correlated wave function theory⁹⁵ (see also Fig. 9), and that we expect the friction coefficients that would be computed at the barrier geometries for O₂ + Al(111) to be small, as the metal electron densities should be small at these geometries, which occur far from the surface.^41,95

Fig. 9 Sticking probabilities computed with the GGA PBE and RPBE DFs,^51,52 the meta-GGA MS-RPBEl DF,⁴¹ the screened hybrid HSE03-1/3X DF,⁴¹ and with embedded correlated wave function theory⁹⁵ are compared with experimental results, for sticking of O₂ on Al(111).⁹⁶ Figure taken from ref. 41 (https://pubs.acs.org/doi/10.1021/acs.jpclett.0c02452). Further permission requests to be directed to the ACS.

Recently two papers^28,58 appeared in which the random phase approximation^99–101 in a DFT framework employing the adiabatic-connection fluctuation-dissipation theorem (ACFDT-RPA,¹⁰² henceforth simply RPA), which may be considered a fifth-rung density functional,²⁰ was used in periodic calculations. The results look very promising: agreement to within chemical accuracy was obtained for H₂ + Cu(111)^28,58 and for H₂ + Al(110)⁵⁸ (see also Fig. 8). Here one has to keep in mind that the method has only been tested on barriers for two dissociative-chemisorption-on-metals reactions thus far. For the BH76 database of gas phase reaction barriers the RPA (MAE = 2.3 kcal mol⁻¹ (ref. 103)) was outperformed¹³ by diffusion Monte-Carlo (MAE = 1.2 kcal mol⁻¹ (ref. 104, and 105)). The reliable application of RPA to energy differences like barrier heights relies on quite substantial error cancellation,²⁰ but this can work quite well in applications to solids.^20,106 The RPA may be viewed as a double hybrid functional and recently another double hybrid functional (XYG3)¹⁰⁷ was tested on H₂ + Cu(111) with very promising results¹⁰⁸ using a cluster-based approach, as will be discussed further below.

Before moving on to electronic structure methods other than DFT it is good to sketch a fundamental problem that exists in DFT for describing reactions of molecules on metal surfaces. Benchmark studies show that for gas phase reactions density functionals exist that exhibit an accuracy for barrier heights of about 2 kcal mol (ref. 23, and 24) To our knowledge, the best DF developed so far²³ is a range separated hybrid DF¹⁰⁹ in which the percentage of exact exchange is maximum at long range (for large distances between the points r₃ and r₄ in Fig. 10). However, for molecule–metal surface systems there is a principal problem: metals are best described with screened hybrid DFs,⁹³ with the amount of exact exchange becoming zero at long range (for large distances between r₁ and r₂ in Fig. 10). In contrast, in the gas phase the fraction of exact exchange would ideally be 1 at long range to correctly describe the required 1/R behavior of the exchange–correlation potential in the gas phase.^110,111 A DF incorporating the correct long-range behavior in the gas and the metal would move us a step closer to a universally accurate DF. A possible route to meeting this challenge will be sketched in Section 3.

Fig. 10 The exchange conundrum of DFT in systems containing molecules and metals. The points r₁ and r₂ are located in the metal, and the points r₃ and r₄ in the gas phase. For additional explanation, see the text.

Diffusion Monte-Carlo (DMC) is a stochastic first principles method.¹¹² It has been applied to DC on metals^{27,37,38,57,113} and to water addition to CO on Pt(111).¹¹⁴ For cases where the minimum barrier to DC on a metal could be compared to experiment (H₂ + Cu(111)²⁷ and Al(110)³⁷) results close to chemical accuracy were achieved (errors of about 1.5 kcal mol⁻¹). The barrier computed for H₂ dissociation above a bridge site on Pt(111) agreed with a semi-empirical SRP-DFT value to 0.8 kcal mol;¹¹³ this comparison is less direct because the bridge site is not the minimum barrier geometry for this system. Agreement with the activation barrier for water addition to CO on Pt(111) (better than 1 kJ mol⁻¹)¹¹⁴ was even more excellent, but perhaps fortuitously so. Comparison of DMC to DFT calculations using standard functionals (with only semi-local exchange) showed that, while generally inaccurate at predicting absolute barrier heights, such functionals are remarkably accurate at predicting how the barrier height varies with the geometry of the barrier in reduced dimensionality.³⁸ This observation helps with explaining why, with tuning of just one parameter in a functional with semi-local exchange, SRP-DFT is usually capable of accurately describing both the threshold and the slope of measured sticking curves.^13,37

For the calculation of gas phase reaction barrier heights, the CCSD(T)¹¹⁵ method represents the gold standard. Calculations on physisorption of small molecules on semi-conductor¹¹⁶ and insulator^117,118 surfaces have been done with the CCSD(T) method, using periodic boundary conditions^117,118 and embedded cluster formalisms.^116,118 These calculations are quite accurate (sub-chemical accuracy). Periodic CCSD(T) calculations of the barrier height for DC of H₂ on Si(100) were consistent with the experimental lower bound of E_b,¹¹⁹ but to our knowledge it has not yet been possible to perform CCSD(T) calculations on molecules interacting with metal surfaces. Zhang and Grüneis have recently reviewed and presented an outlook for the application of coupled cluster methods in the field of materials science, including a discussion of surface chemistry.¹²⁰

CCSD(T) but also other ab initio many-electron wave function methods are currently too computationally expensive for periodic calculations on DC on metals. However, they may be applied in an embedded cluster fashion (“embedded correlated wave function theory”, ECW), using density functional embedding.¹²¹ While this leads to a formally exact theory,¹²¹ in practice a limiting factor is the size of the cluster that can be treated with ab initio theory, as the convergence of the results depends critically on the size of the cluster that can be handled.^26,28 Nonetheless accurate results can be obtained with this approach. Dynamics calculations based on an embedded CASPT2 (ref. 122 and 123) PES were in semi-quantitative agreement with experiments on sticking of O₂ on Al(111),⁹⁵ although agreement to within chemical accuracy was not yet achieved (see also Fig. 9). However, embedded CASPT2 calculations failed for H₂ + Cu(111), which at the time was not recognized because the comparison to experiment was not made in an appropriate way:¹²⁴ it was based on an erroneous fit of the activation energy as a function of temperature¹²⁵ instead of a direct comparison to the semi-empirical barrier height.³³ On the other hand embedded NEVPT2 (ref. 126) calculations agreed with the semi-empirical barrier height for H₂ + Cu(111) to within chemical accuracy,²⁸ and there is reason to believe that NEVPT2 should be more robust than CASPT2.¹²⁴ In the same way as used for correlated wave function theory density functional embedding can also be used with RPA for the cluster, which provided semi-quantitative agreement with the semi-empirical barrier height for H₂ + Cu(111) at far less computational expense than more accurate periodic RPA calculations.²⁸ A problem with validation of results may be that calculations in which the molecule is centered on different clusters for different adsorption geometries in different ways might be needed, while using these calculations to construct a PES in a consistent manner may then not be obvious.

Another approach in which clusters are employed to allow the use of correlated wave function methods is the ONIOM approach.¹²⁷ In this approach the ab initio result for the molecule interacting with the metal cluster is corrected with the difference of a periodic and a cluster DFT calculation for the same adsorption geometry. Sautet and co-workers¹²⁸ applied this approach to H₂ dissociation on Cu(100) at a barrier geometry in reduced dimensionality (this was not the full-dimensional minimum barrier geometry). Good agreement with a semi-empirical SRP-DFT result was obtained with the use of MRCI-Q (i.e., multi-reference CI¹²⁹ with Davidson correction¹³⁰) but not with a CCSD(T)¹¹⁵ treatment of the cluster, which seems suspect. Araujo et al.¹³¹ have used the ONIOM approach employing a global hybrid (i.e., M06) density functional⁹¹ for the cluster, while using a dispersion corrected DFT functional (PBE-D3 (ref. 80, 132 and 133)) to correct for finite-size effects. One might argue that in such a set-up a hybrid functional is used that is screened in a primitive way, with the screening length determined by the cluster size. Results¹³¹ were obtained that were in excellent agreement with accurate semi-empirical reference values of minimum barrier heights for five systems contained in the SBH10 (ref. 31) and SBH17 (ref. 32) databases. However, the reference values were, for 4 of the 5 systems studied, classical barrier heights obtained with the surface held static at the geometry appropriate for the vacuum–solid interface. Because Araujo et al.¹³¹ computed activation energies (i.e., zero-point vibration energy corrected barrier heights, allowing the surface atoms to relax in the calculations on the transition state), the agreement that was achieved is rather meaningless.³⁷ Xu and coworkers¹⁰⁸ have made this important distinction when benchmarking the performance of the doubly hybrid functional XYG3 (ref. 107) for H₂ dissociation on Cu(111) using an ONIOM-based approach including a systematic convergence test with respect to cluster size. It is very encouraging that they have obtained remarkably good agreement (to within chemical accuracy) with the semi-empirical SRP-DFT result.¹⁰⁸

In summary, the state-of-the-art for calculating classical barrier heights for DC on metals is as follows. DFT with functionals using semi-local exchange yields a mean absolute error ≥2.4 kcal mol⁻¹ for easy-to-model systems. These functionals tend to be accurate for the variation of the barrier height with system geometry. As a result, tuning such functionals to obtain agreement with experimental sticking curves (SRP-DFT) has facilitated the construction of a database (SBH17) with chemically accurate reference values of classical barriers for easy-to-model systems. However, functionals with semi-local exchange systematically underestimate barrier heights for difficult-to-model systems. Calculations on O₂ + Al(111) offer hope that screened hybrid functionals are more accurate for these systems. Limited evidence for two H₂-metal systems (H₂ + Cu(111) and Al(110)) offers hope that this increased accuracy holds for molecule–metal surface systems in general. Evidence for these two systems suggests that the most accurate approaches currently available are periodic DMC and periodic RPA. Correlated wave function theory with density functional embedding is accurate in principle but problems with the convergence of the size of the cluster treated with the high-level method may stand in the way of achieving chemical accuracy. This cluster size convergence problem does not occur for calculations on H₂ + Cu(111) with the doubly hybrid functional XYG3 based on the ONIOM approach.¹⁰⁸ Periodic calculations with these functionals and without embedding can hopefully confirm this very encouraging trend.

2.2 State-of-the-art in describing electronically non-adiabatic effects in molecules scattering from metal surfaces

In periodic solids the electrons form a continuum of states (bands). In metals, the highest band with electrons in it, the conduction band, is only partly filled. With unoccupied states readily available at practically zero excitation energy, it follows that collisions with molecules can (or actually, will^134,135) at least cause transitions in which electrons are excited from states below the metal's Fermi level to states above the Fermi level. This process is called electron–hole pair (ehp) excitation, and it breaks the Born–Oppenheimer approximation. Globally speaking, two cases can be distinguished.¹³⁶ In the first case scattering takes place without changes in the electronic state of the molecule scattered back to the gas phase, or, in case of reaction, with the electronic state of the system staying close to the electronic ground state of the full system. In the second case, the molecule scattered back to the gas phase has changed its electronic state, or it has done so temporarily when it interacted with the metal, or immediately after reaction the electronic state of the system is not close to that of the ground state of the full system. Loosely speaking the first case might be thought of as a case of weak coupling¹³⁶ between the nuclear and electronic degrees of freedom, while in the second case there is strong coupling¹³⁶ between the two. Below we will focus on the first case but we will also say something about the second case. We will restrict ourselves to methods that have been applied to actual chemical/physical systems, although excellent work (e.g. ref. 137–139) has been presented concerning methods that have been applied to limiting case model systems, which models in future may also be applied to real life systems.

Systems with weak coupling may be dealt with using electronic friction (EF) methods.^66,140 In the words of Dou and Subotnik,¹⁴¹ “EF represents the first order correction to the Born–Oppenheimer approximation in the presence of a manifold of fast electronic states”, i.e., those of the perturbed metal. As a result, atoms moving close to a metal surface “experience a drag in the presence of a manifold of electronic states”.¹⁴¹ The EF method may be seen as an extreme type mean field Ehrenfest method, in which the electronic ground state potential replaces the potential averaged over (multiple) ehp states. In MDEF, the energy dissipation to the electrons directly affects the nuclear motion. Time-dependent perturbation theory can be used to derive electronic friction theory, which can be cast in a generalized Langevin form for the nuclear equations of motion.⁶⁶ Friction coefficients mediate the energy flow to the electrons through friction forces, and through fluctuation forces that depend on temperature. EF theory invokes the Markov approximation, thereby assuming short (i.e., zero) decoherence times.^136,142

Basically two types of EF methods exist. In orbital dependent friction (ODF)^66,143–146 the elements of the EF tensor (i.e., the friction coefficients) are written as

The friction coefficients are elements of a N_d × N_d friction tensor. Here N_d is the number of atoms subjected to electronic friction times 3 for the number of Cartesian directions, so for a diatomic molecule subject to friction N_d = 6. In eqn (2)i and j refer to atoms i and j, α and β to Cartesian directions, and g_kab^iα(R) is an electron–phonon coupling matrix element. For the motion of adsorbate atom i along direction α g_kab^iα(R) describes the resulting non-adiabatic coupling between two electronic states with band indices a and b at wave vector k. These two electronic states are states of the whole system (the molecule plus the metal). The electronic structure of the molecule as well as that of the metal is therefore taken into account in ODF theory,^66,143–146 and this is an advantage of the ODF method. Finally, ε_ka is the energy of the electron with band index a at wave vector k, and ε_Fis the energy of the Fermi level.

ODF coefficients generally yield anisotropic friction tensors that may depend strongly on the molecular coordinates. In performing ODF calculations one has to take special care that the computed friction coefficients are well converged with respect to a number of parameters in the underlying electronic structure calculations.¹⁴⁷ For computationally convenient evaluation in dynamics, friction tensors can be fitted using symmetry adapted neural network techniques.¹⁴⁸

An advantage of the ODF method already mentioned above is that it takes into account the electronic structure of the molecule and the metal. However, the ODF method also has several disadvantages. ODF is based on linear response theory, which means that non-adiabatic perturbations are taken into account only up to first order, which further emphasizes the character of the method as a “weak coupling” approximation.^66,149 Periodic DFT calculations of ODF tensors (i.e., the computation of the friction coefficients in eqn (2)) for an isolated molecule interacting with a surface are typically based on inter-band transitions only, combined with a broadening of the δ-functions present in eqn (2) during Brillouin zone integration. This leads to a non-physical dependence of the resulting friction coefficients on broadening parameters,^69,150–152 which represents a conceptual shortcoming. In practice the problem has been handled pragmatically by performing calculations with varying broadening parameters. Often one can show that the results are rather insensitive to these widths, within a range of widths that includes the broadening parameter used in the actual calculations.^146,153 Nonetheless this is not a wholly satisfactory state-of-affairs. ODF has also been shown to yield unphysically large friction in the region of the doublet-singlet spin transition for H-atoms incident on a non-magnetic metal surface (Cu(111)).¹⁵⁴ Because techniques were lacking for the efficient evaluation of eqn (2), until recently ODF coefficients were only computed for and used in low-dimensional scattering calculations^155–158 However, recently 6D friction tensors have become available for and have been used in dynamics simulations on DC of diatomic molecules on metals (e.g., H₂ + Ag(111)^152,153 and Cu(111),¹⁴⁶ and N₂ + Ru(0001)⁶⁹).

The local density friction approximation (LDFA) is the second state-of-the-art EF technique currently in use. The LDFA maps the friction problem onto single atoms embedded as impurity in a homogeneous electron gas (HEG) via the electron density of the pristine, unperturbed surface at the points where each atom in a molecule resides.¹⁵⁹ The HEG, or jellium, is a simple model for bulk metals with nearly free electrons in which the nuclei of the metal atoms have been smeared out as a constant background. The scattering of jellium electrons at an atom embedded therein has been calculated at the level of the local density approximation (LDA)¹⁶⁰ and more recently also at the generalized gradient approximation¹⁶¹ to DFT by self-consistently solving the Kohn–Sham equations, using spherical symmetry and atomic units

for this atom-in-jellium (AIJ) system.¹⁶⁰v^AIJ(r) is the effective scattering potential by which the atom perturbs the jellium background, including metallic screening. It is defined by the element type of the atom and the jellium electron density. Electrons at the Fermi level with momentum and angular momentum l, described by radial wave functions ψ_l(r) scatter from this potential, which results in the well-known phase shift expression for asymptotically large distances r = R¹⁶¹In eqn (4), the j_l and n_l are spherical Bessel and Neumann functions, respectively. The δ_l can be obtained from the literature for a large amount of elements and jellium densities.^160,161 In the LDFA the EF coefficients are then obtained from the phase shifts using

LDFA friction coefficients rely on the independent atom approximation (IAA): they are isotropically dependent on the Cartesian coordinates of each atom in the molecule.¹⁵⁹ Thus, a disadvantage of the LDFA method is that the effect of the full ES of the molecule is not taken into account.¹⁶²

Eqn (2) and (5) can be used to compute friction forces that are linearly proportional to the atom's velocity and the friction coefficients defined in these equations. The LDFA rests on a firm theoretical basis when applied to atoms scattering from metals.¹⁶³ The approach based on eqn (5) and the friction force derived from it have been shown to correspond to the exact low velocity limit of an atom moving through jellium as obtained with time-dependent density functional theory.¹⁶⁴ The LDFA has been applied to many DC systems^{69,146,152,153,165,166} and much of the work prior to 2017 has been reviewed by Alducin et al.¹⁹

As discussed also in ref. 13, calculations with the LDFA by itself suggest that ehp excitation can usually be neglected when computing DC probabilities, with energy shifts between adiabatic and non-adiabatic sticking curves usually smaller than 1 kcal mol⁻¹. Comparisons now available for H₂-metal systems (H₂ + Cu(111)¹⁴⁶ and Ag(111)^152,153) also suggest that for these systems it hardly matters whether ODF or the LDFA is used to compute EF coefficients. The only system for which a large difference has been found between sticking probabilities computed with the LDFA and ODF is the already mentioned N₂ + Ru(0001) system (see Fig. 5).⁶⁹ Much of the above observations can be explained from the size of the friction coefficients for the systems discussed, the differences between LDFA and ODF values (Fig. 11), and the regime of velocities relevant to the sticking. For instance, the ODF friction coefficients for N₂ + Ru(0001) are roughly an order of magnitude larger than those for H₂ + Cu(111) (Fig. 11). Also, for N₂ + Ru(0001) the friction coefficients for motion towards the surface as well as for the molecular vibration are much larger with ODF than in the LDFA. The same is not true for H₂ + Cu(111). Even though N₂ is much heavier than H₂, their velocities do not differ much in the regimes relevant to reaction (up to about 0.5 eV for H₂ + Cu(111),¹⁴⁶ and 4 eV for N₂ + Ru(0001)⁶⁹). The much larger ODF friction coefficients for both degrees of freedom for the latter system then make the difference, explaining both why ODF friction has a large effect, and why this effect is much larger than obtained with the LDFA. These results also suggest caution in interpreting the older LDFA results: the N₂ + Ru(0001) results of Fig. 5 convincingly show that a small effect of LDFA friction on sticking is no guarantee whatsoever for a small effect of ODF on sticking (see Fig. 5). This latter observation makes research on the accuracy of these two EF theories, and efforts aimed at developing a better theory all the more relevant.

Fig. 11 Diagonal elements of the friction tensor, η_qq(in m eV ps Å⁻²), are shown as a function of q (d is the bond distance of the molecule, Z the distance of the molecule's center of mass to the surface) for the LDFA and the ODF approximations, as computed for H₂ + Cu(111)¹⁴⁶ and N₂ + Ru(0001).⁶⁹ Data taken from ref. 146 and 69. This figure has been reproduced from ref. 13 with permission from the Royal Society of Chemistry, copyright 2021.

The IESH method has originally been developed for the strong coupling regime, where EF approaches are not expected to work.¹⁶⁷ Perhaps the most prominent example is multi-quantum vibrational relaxation of initially highly vibrationally excited NO scattering from Au(111) (see Fig. 3d). The IESH method uses a discretized version of the Newns–Anderson model, in which electrons do not interact.^67,68,167 In these one-electron states the modeled electrons can be present in states below and above the Fermi level in the conduction band of the metal, and in the affinity level, i.e., the extra orbital of the molecule to which a metal electron can move to make a state in which the anion of the molecule interacts with the metal (a π* orbital in NO). The method can thus be used to simulate a system in which (partial) electron transfer from the surface to the molecule occurs. The N_e electrons included in the model can be in the lowest M one-electron states; then the system is in the ground state. Alternatively, one or more electrons can be excited to higher-lying one-electron states. The motion of the nuclei is treated with the QCT method subject to the forces related to the PES for the electronic state the system is in at any given time. The time-dependent Schrödinger equation (TDSE) for the electrons is integrated simultaneously. The effect of the NA couplings (between the neutral molecule–surface state and the anion–surface state, and through ehp excitation) is taken into account using the fewest switches surface hopping method.¹⁶⁸

Important inputs to the IESH calculations are the PESs for the diabatic neutral molecule–metal and anion–metal states, and the coupling potential. In the approach used mostly so far⁶⁵ they have been obtained from the adiabatic ground state PES, the Bader charge of the molecule, and the change of the electronic energy of the system upon application of a small electric field. Meng and Jiang¹⁶⁹ have recently described a pragmatic approach in which two diabatic PESs are computed with constrained DFT (CDFT^170,171), constraining both the surface and the molecule to be neutral in one case, and putting an extra electron on the impinging molecule in the other case. Hereby they used that with CDFT the number of electrons on the molecule can be set equal to the number appropriate for either the neutral molecule or the anion, thereby obtaining the appropriate corresponding diabatic energy.^170,171 CDFT can also be used to derive the coupling potential,¹⁷¹ but for this Meng and Jiang instead followed the original approach of Tully and co-workers.⁶⁵ Diabatic PESs and couplings were analyzed to arrive at qualitative predictions regarding vibrational relaxation of NO and CO scattered from Au(111) and Ag(111).¹⁶⁹ The revPBE DF¹⁷² was used in the calculations on NO + Au(111) and Ag(111), and the vdW-DF1 DF⁷⁵ (featuring revPBE exchange but replacing PBE with vdW-DF1 correlation) for CO + Au(111) and Ag(111).

Very recently Jiang and co-workers¹⁷³ have used their CDFT approach to construct high-dimensional neural network PESs for CO + Au(111). They have used their approach to study vibrationally inelastic scattering of initially highly vibrationally excited CO (in v = 17) from Au(111). As can be seen in Fig. 12 their IESH calculations¹⁷³ built on CDFT potentials yields results in quite good agreement with the original experiments.¹⁷⁴ However, the improvement over dynamics calculations using the same PES but making the Born–Oppenheimer approximation (“BOMD”) is not unequivocal. For example, while the BOMD results overestimate the survival probabilities of CO(v = 17), the IESH results tend to underestimate these quantities. The improvement achieved with IESH for vibrational de-excitation from v = 2 to v = 1 is however unambiguous (Fig. 4a of ref. 173, results not shown here). Nonetheless one should keep in mind that reproducing similar results for NO scattering from Au(111) (as in Fig. 3) should put higher demands on the approach to compute PESs, as the charge transfer energy of NO + Au(111) (5.3 eV) is far lower than that of CO + Au(111) (6.9 eV).

Fig. 12 Measured¹⁷⁴ relative probabilities for scattering of CO from Au(111) in its initial vibrational state v_i = 17 to the final vibrational state v_f (on the x-axis) are compared with calculations¹⁷³ using the IESH method. For details see the text. Taken from ref. 173. Reprinted figure with permission from [G. Meng, J. Gardner, N. Hertl, W. Dou, R. J. Maurer and B. Jiang, Physical Review Letters, 133, 036203, 2024]. Copyright (2024) by the American Physical Society.

2.3 State-of-the-art tools for dynamics

Dynamics calculations can be used to calculate observables such as sticking probabilities and probabilities for vibrationally and/or rotationally inelastic and/or diffractive scattering.¹³ The hole model³⁹ teaches us that agreement of computed with measured sticking probabilities indicates that the barrier heights obtained with the electronic structure method used are accurate.^13,40 This is the case if direct dynamics or an accurately fitted PES is used, and if an accurate dynamical model and method is selected to compute the sticking probabilities. Good agreement of dynamics calculations with experiment for more detailed dynamical observables, such as rotationally inelastic scattering probabilities or diffraction probabilities, indicates a very high quality of the underlying electronic structure method and possibly the method used to fit the PES, as such computed observables are highly sensitive to details of the PES used.¹⁷⁵

The fitting of a PES can be avoided if one can afford the computational resources to run direct dynamics calculations (ab initio molecular dynamics (AIMD)^176,177 or density functional molecular dynamics (DFMD)³⁴), which uses (quasi-)classical trajectories. This approach will work for the calculation of probabilities that are not too small (currently, >10⁻³), as the inherent statistical error in the computed probability p is equal to , where N_t is the number of trajectories run. Otherwise, it is best to first compute electronic structure data and fit a PES to these data. Fitting PESs can now be done accurately with the corrugation reducing procedure¹⁷⁸ for diatomic molecules interacting with static surfaces, using permutationally invariant polynomials combined with neural network methods for diatomic and polyatomic molecules interacting with static surfaces,^179,180 or using a high-dimensional neural network method^70,181,182 for molecules interacting with mobile surfaces, in which the surface atoms are allowed to move. In general, fitting an accurate PES to electronic structure data does not represent a bottleneck to the accurate calculation of observables. However, efficiency may be an issue if there is high computational expense associated with either the electronic structure calculations (cost per data point) or the dynamics calculations (number of trajectories needed).

The choice of the dynamical model requires care. As a rule of thumb, for accurate sticking probabilities surface atom motion needs to be modeled if the impinging molecule is heavier than D₂, or the surface temperature is high (T_s ≫ 300 K).¹³ In another rule of thumb, it is advisable to attempt to model electronically non-adiabatic effects if E(CT) is smaller than or close to 7 eV (ref. 69) and/or the impinging molecule is in a highly excited vibrational state.^15,183,184 The role of both dissipative mechanisms can now be studied with MDEF calculations if a high-dimensional neural network PES is available,^44,69,185 or with ab initio molecular dynamics with electronic friction (AIMDEF).^186,187 As noted above, the accuracy likely depends on the type of EF implemented, but it is not yet clear which approach is most accurate. Strong non-adiabatic effects can be modeled with the independent electron surface hopping (IESH) method⁶⁸ using classical dynamics and, essentially, the fewest switches surface hopping algorithm,¹⁶⁸ and this can now be done while also including surface motion accurately.¹⁷³

The quasi-classical trajectory (QCT) method is, perhaps surprisingly, usually very accurate for modeling activated DC of diatomic molecules as measured in supersonic molecular beam experiments, even for H₂.¹⁶⁵ This is because¹⁸⁸ these experiments imply a high degree of averaging over rovibrational states, and, especially for highly activated dissociation, in experiments at low incidence energies sticking often takes place through the (classically allowed) reaction of vibrationally excited molecules.^189,190 Also for this reason tunneling usually does not play a very large role.¹⁹¹ However, in DC of polyatomic molecules intramolecular vibrational relaxation may cause problems for the QCT method, especially for molecules incident in vibrationally excited states and/or at high nozzle temperatures.¹⁹² These problems may be avoided by e.g. restricting simulations to cases where a molecule containing H-atoms is partially deuterated with only the one remaining XH-stretch vibration pre-excited, or by keeping the nozzle-temperature low.^35,192,193 In the simulation of state-to-state scattering experiments special attention has to be paid to the algorithm to assign final states.^194–197 With the use of particular Gaussian binning methods for assigning final states and a method called adiabatic correction it is also possible to obtain QCT sticking probabilities in good agreement with quantum dynamics results for non-activated sticking, as shown by Bonnet and co-workers for H₂ + Pd(111).¹⁹⁶

Quantum dynamical simulations of DC of diatomic molecules have usually been performed with the time-dependent wave packet method,^198–201 the application of which is fairly routine by now if the molecule contains at least one H-atom. However, TDWP calculations on reaction of heavy diatomic molecules on static surfaces can still be difficult to obtain converged results for, and are rather rare.²⁰² In the largest DC problem to which the TDWP method was applied while simulating motion in all molecular degrees of freedom without approximations, H₂O reacts on Cu(111)²⁰³ and Ni(100).²⁰⁴ Calculations with the TDWP method on DC of polyatomic molecules often invoke reduced dimensionality approximations.^54,205,206 In TDWP calculations usually the static surface approximation is invoked, but reliable approximations have been developed to take into account the effect of surface atom motion or surface temperature on sticking on metal surfaces in an a posteriori fashion.^207–209

The DC of polyatomic molecules, like CH₄, H₂O, and CO₂, on metals has often been studied with a reaction path Hamiltonian (RPH) method.^210,211 While approximate, a large advantage of this method is that the effect of all the molecular vibrations can be modeled for such small-sized polyatomic molecules. While approximations are usually needed to the molecular rotations and translational motion parallel to the surface, guidelines on successful strategies are available.^{205,212–215} An already mentioned a posteriori method^207,208 can be used to accurately^216,217 treat the effect of surface temperature. Under conditions in which quasi-classical mechanics should be accurate RPH calculations agreed well with full-dimensional DFMD calculations, where all calculations were based on the same density functional.^191,214 However, it has been suggested that under conditions at which the system may swerve from the minimum energy reaction path (e.g. at high incidence energies) the RPH method may become less accurate.¹⁹²

A fairly new approach to incorporating quantum effects in molecular dynamics simulations is non-equilibrium ring polymer molecular dynamics (NE-RPMD).²¹⁸ Strictly speaking RPMD is only valid under specific conditions not found in real systems,²¹⁹ and just like the QCT method it therefore needs to be tested on systems to establish under which conditions it is valid. Comparisons to quantum dynamics results for e.g. H₂ + Cu(111),²²⁰ D₂O + Ni(111),²²⁰ and H₂ + Co(0001)²²¹ show that NE-RPMD performs much better than the QCT method in the tunneling regime, probably also because NE-RPMD is capable of avoiding artificial energy leakage. However, the method did not perform as well as QCT for H₂ + Cu(111)²²⁰ at intermediate collision energies, and slightly overestimated reaction in H₂ + Co(0001) at high incidence energies.²²¹ A comparison to quantum dynamics results for H₂ + Pd(111)²²¹ shows that NE-RPMD may not work for non-activated dissociation if the dynamics is affected by quantum resonances. An advantage of NE-RPMD is that its use avoids artificial intramolecular vibrational redistribution (IVR), making it a much better method than QCT for describing DC of CH₄ on Pt(111).¹⁹² The ordinary RPMD method applicable in the canonical regime performed rather well at reproducing accurately measured experimental rate constants for DC of H₂ on Pt(111).²²² This finding was significant because the calculations²²² were based on a SRP density functional fitted to supersonic molecular beam experiments,⁸⁶ which should therefore yield a chemically accurate description of the H₂–Pt(111) interaction.¹³ NE-RPMD has also been used to study scattering of H-atoms from graphene,²²³ hydrogen spillover from a Pt atom embedded in Cu(111),²²⁴ and NO desorption from Pd(111) under thermal conditions.²²⁵

3. The way forward

It is clear that the largest challenges exist for systems characterized by E(CT) < 7 eV. Therefore, these are the systems to tackle. To be practical, it is good to pick systems that have already been studied experimentally, because especially in the beginning it will be necessary to validate computed results through comparisons with experiments. Systems with E(CT) < 7 eV for which accurate results for activated sticking are available from supersonic molecular beam experiments include O₂ + Al(111),^96,226 Ag(110),^227–229 Cu(111),^230,231 and Cu(100),^232–234 HCl + Au(111)⁴³ and Cu(111),²³⁵ D₂O + Ni(111),⁵⁹ CO₂ + Cu(110),²³⁶ and NO + Cu(111).^237,238 There are thus plenty of systems to study. Experimental results for state-to-state scattering are also available for HCl + Au(111)^239,240 and for NO + Cu(111).^237,238 As will be discussed below the availability of experiments on sticking and state-to-state scattering has special advantages in the situation where uncertainties exist concerning the electronic structure method as well as the method to treat non-adiabatic energy dissipation.

3.1 Electronic structure

As mentioned before the greatest challenge is to come up with an electronic structure method that is accurate and ultimately based on first principles as well as efficient at computing PESs. For this we suggest to use parameterized density functionals (similar to eqn (1)) as before, but now to base the choice of functional on a few calculations with an accurate, first principles-based (FPB) electronic structure method for judiciously chosen (barrier) geometries. Thereby FPB-DFs are still specific to a particular system, but they do provide a steppingstone towards broader testing of the underlying first principles electronic structure methods. An approach that uses the diffusion Monte-Carlo method to generate the first principles results (quantum Monte-Carlo based DFT, QMC-DFT) was recently tested successfully on the DC of H₂ on Al(110).³⁷ The proof-of-concept calculation provided a sticking curve that was displaced from the measured sticking curve²⁴¹ along the incidence–energy axis by only 1.4 kcal mol⁻¹, suggesting the barriers in the PES to be accurate to within this number (see Fig. 13). For this system E(CT) = 7.2 eV, and a semi-local exchange functional could be used that was a weighted average of the exchange parts of the PBE⁸⁰ and RPBE⁴⁸ GGA DFs, while vdW-DF2 (ref. 76) correlation was employed. The QMC-DFT approach was called a best-of-both-worlds approach³⁷ because it combined the accuracy of DMC (for the few calculations “nailing the PES” to first principles results) with the efficiency of DFT for computing a whole PES. Just like the SRP-DFT approach, QMC-DFT works because standard DFT is already quite good at describing the variation of the barrier height with system geometry,^37,38 so the only thing that needs to be done is to fit the parameter in the DF to the first principles result for the minimum barrier height.

Fig. 13 Computed³⁷ and measured^241,242S₀ for H₂ + Al(110). S₀ computed with the NBOMS model, and as computed with the NBOMS model but also corrected for quantum dynamical effects, are compared to experimental values. The blue horizontal lines and numbers (kcal mol⁻¹) indicate the energy distance between the measured S₀ and the S₀ computed with the NBOMS model corrected for nuclear quantum effects. Figure taken from ref. 37 (https://doi.org/10.1021/acs.jpclett.3c02972). Further permission requests to be directed to the ACS.

The QMC-DFT approach is yet to be tested on systems with values of E(CT) less than 7 eV. For reasons discussed above, we anticipate that the density functional of eqn (1) will have to be modified to include exact exchange (EXX). Based on the research on O₂ + Al(111)⁴¹ discussed above such a DF should combine screened hybrid exchange with van der Waals correlation. To reiterate, measured S₀ for O₂ + Al(111) could be reproduced semi-quantitatively with the HSE03-1/3X DF,⁴¹ which is the screened hybrid HSE03 DF^93,94 with the maximal fraction of EXX, a, increased to 1/3 (see Fig. 9). Increasing a further and adding van der Waals correlation should broaden the computed S₀ curve,⁴¹ increasing the agreement with the experiments⁴¹ and a similar approach is likely to work for other systems with E(CT) < 7 eV. An even better approach will probably be to use one of two recent screened hybrid van der Waals DFs, called vdW-DF-ahcx²⁴³ and vdW-DF2-ahbr,²⁴⁴ which were designed to not only give a good description of metals but also to make the exchange and correlation parts of the functional compatible. The parameter a (and/or the screening range parameter^243,244) can be tuned to achieve agreement with QMC. If these DFs^243,244 do not work for a specific system, this may not present a problem: there are sufficient other conceivable combinations of screened hybrid^111,245 and van der Waals DFs^77,246 to try. We also believe that the balance that can be struck between the need to correct for the self-interaction error and the need to describe static correlation by tuning a will not depend much on the geometry of the barriers relevant to the systems we will investigate. Should this not be so then a way to correct for this exists (see below).

We believe that DMC (as primed with a variational Monte-Carlo calculation using a single Kohn–Sham DFT determinant as input^27,37) is, of the first principles methods already tested, likely the most successful. However, the available evidence^28,58,108 suggests that the RPA^99–102 method and the double hybrid functional XYG3 (ref. 107) should also be tried. In this respect a more general name for the approach than QMC-DFT could be first principles-based DFT (FPB-DFT), on the understanding, of course, that essentially all DFs are based on first principles via LDA correlation. Should it be possible to extend CCSD(T)¹¹⁵ from molecules interacting with semi-conductor surfaces¹¹⁹ to molecule–metal surface interactions, then naturally CCSD(T) can also be applied within an FPB-DFT approach. The same goes for the ECW method of Carter and co-workers,^28,121 which would be especially worthwhile if the embedded cluster can be made large enough for convergence in combination with a highly accurate correlated wave function method. In all cases the principle would be the same: the accuracy would come from fitting a mixing coefficient in a density functional of an appropriate form to accurate first principles results for well-chosen geometries, and the resulting density functional would then allow the efficient computation of a full PES.

Coming back to DMC, one may attempt to improve DMC over the single-determinant version tested so far. The main sources of error in a DMC calculation are the fixed-node approximation (the quality of the DMC solution depends somewhat on the quality of the trial wave function through the nodes it inherits from it) and the locality approximation, which arises from the need to use pseudo-potentials to keep the costs of the calculation at bay. Attempts to achieve chemical accuracy (errors <1 kcal mol⁻¹, chemical accuracy) would probably need to reduce the fixed node error using a multi-determinant approach or approaches equivalent to it and based on geminals.^247,248 A multi-determinant approach one could try would use charge constrained DFT^170,171 to compute diabatic DFT wave functions for the neutral system and the system with one electron transferred to the molecule,¹⁶⁹ using that the two Slater determinants used need not be orthogonal.²⁴⁹ For systems containing 5d metals like HCl + Au(111) it may be necessary to use a new type of more accurate pseudo-potential, but procedures to obtain these^250,251 or more general procedures to reduce the locality error²⁵² are available. It may also be possible to reduce the statistical error in the E_b computed with DMC using approaches aimed at reducing time step error.²⁵³ In other words, there are systematic ways of improving the quality of the DMC calculations used to pin the barrier height in QMC-DFT.

In Section 2.1 we have already sketched a fundamental problem existing with the treatment of molecule–metal surface reactions using DFT, also referring to Fig. 10. It would be nice to have a generally applicable DF (i.e., not one specific to only one or a few particular systems, like a SRP-DF¹³ or a FPB-DF) available with a guaranteed accuracy of ≈ 2 kcal mol⁻¹ for barrier heights for DC of molecules on metal surfaces. As sketched above, for our system of interest in the gas phase the DF would ideally behave as a range separated hybrid DF¹⁰⁹ in which the percentage of exact exchange is maximum at long range (for large distances between r₃ and r₄ in Fig. 10). However, in the metal the DF should behave as a screened hybrid DF,⁹³ with the amount of exact exchange becoming zero at long range (for large distances between r₁ and r₂ in Fig. 10). One way to accomplish this would be to develop a “true made-simple hybrid DF”.^254,255 To understand this class of DFs we draw an analogy with “made-simple” meta-GGA DFs.^79,88 In these DFs, in each point of space an optimal mixture of two GGA exchange DFs is determined from an inhomogeneity parameter α that depends on the local kinetic energy density τ.⁸⁸ The value of α is used to determine whether the τ computed at a point in space corresponds to metallic or covalent bonding.^88,256 Extension of this idea to hybrid DFs implies that the percentage of exact exchange assigned to each pair of points in the double integral over space should depend on the values of α at both points, leading to what we call a true made-simple hybrid DF (which differs from the made-simple hybrid DF introduced in ref. 88). It should also be possible to use the values of other parameters²⁵⁶ at these points.²⁵⁷

To arrive at a true-made simple hybrid functional, in evaluating EXX in the integral over the two points r and r′ one would like the contribution to the exchange made by these points to depend on where they are in Fig. 10. For this we first note that the expression for a global hybrid DF reads:

E_xc,gh = aE_x,gh + (1 − a)E_x,sl + E_c (6)

Here, E_x,gh is the EXX DF, E_x,sl is a semi-local exchange DF, and E_c is a correlation DF. In a MS meta-GGA DF,²⁵⁸ the fraction of EXX, a, is set to zero, and the exchange enhancement factor defining E_x,sl is written as

F_x,int(s,α) = F_x,1(p) + f(α)[F_x,0(p) − F_x,1(p)] (7)

In eqn (7), α is a function of the kinetic energy density τ, and α and f(α) are constructed in such a way that f(α) = 0 corresponds to the uniform electron gas and f(α) = 1 to a single-orbital regime.⁸⁸ By choosing F_x,1(p) and F_x,0(p) to be good DFs for metals and molecules respectively, Sun et al. designed a meta-GGA DF that is good for atoms, molecules, surfaces, and solids,⁸⁸ and this design principle has been used to construct meta-GGAs that work well for molecule–metal surface systems with E(CT) > 7 eV.⁷⁹

Sun et al. then went on to define what they called a MS hybrid DF.⁸⁸ However, instead of making a a function of f(α), which might have seemed the logical choice, they merely used f(α) to define the gradient enhancement factor in eqn (7), and then fitted a to a database, this way defining a global hybrid DF, which they called MGGA_MS2h, and which had a = 0.09.⁸⁸ We would not expect this DF to perform well for barrier heights for DC on metals in general, because a DF with this low a value of a would not be expected to work well for systems with E(CT) < 7 eV.

Instead, one might attempt to develop a true MS hybrid DF of the form

In eqn (8), n(r) is the electron density at r, s(r) is the scaled gradient of the density at r, and F_x,sl is the exchange enhancement factor for a suitably chosen semi-local exchange DF. This form is similar to a form used recently by Borlido et al.,²⁵⁵ and bears similarities to the form adopted for the PSTS DF.²⁵⁴ For the non-local mixing function g[α(r),α(r′),|r − r′|], inspired by Borlido et al.²⁵⁵ we suggest to test expressions that are symmetric and separable like

For h(α) one could try an expression that is proportional to the f(α) used in the MS meta-GGAs, which approaches zero in metals, where it will thus appropriately suppress EXX. With the scheme outlined, in the integral over r and r′ the contribution to the exchange made by these points would depend on where they are in Fig. 10, as needed. We sketch just the basic ingredients the method might have here. The true MS hybrid DFs bear a similarity to the local hybrid DFs discussed in ref. 257, which use a local mixing function in the first term on the rhs of eqn (8); one might be able to take inspiration and guidance from this work, for instance, concerning the use of a so-called calibration function one will need, which was omitted by Borlido et al.²⁵⁵ Clearly, work will be needed on g, and on the relationship between g and F_x,sl. As noted by Borlido et al.²⁵⁵ it should also be possible to use off-diagonal mixing in the second term on the rhs of eqn (8), which would however imply the use of a non-local exchange DF instead of the semi-local DF used in this term.²⁵⁹

An alternative approach one might also pursue would be to develop a general purpose DF through machine learning. This might simply be a follow-up version of the DM21 local hybrid DF, for which code and learned network weights are publicly available,²⁶⁰ which should then also be trained on e.g. reaction barrier heights for DC on metals and adsorption of molecules on metals. Such a DF would be a local hybrid DF, as the exchange energy density is already an ingredient of DM21. It might also be possible to develop a machine learned DF in the spirit of the made simple hybrid DF approach by providing the integrand of the first integral in eqn (8) on a two-dimensional grid of points (r, r′) to be used as an additional “feature”²⁶⁰ in the training of the underlying neural network.

3.2 Non-adiabatic effects

As discussed above the ODF and LDFA EF methods both have their specific advantages and disadvantages for describing non-adiabatic effects on reaction. To combine their advantages while avoiding their disadvantages, a new scheme we call scattering potential friction (SPF) could work as follows. As discussed before, the LDFA scheme relates the atom-in-jellium model and the “real surface system” via the electron density in the bare metal. Instead, the SPF scheme would extract an electronic scattering potential from a DFT slab calculation that includes the full molecule–surface interaction, unlike the one obtained self-consistently in the atom-in-jellium model. This Kohn–Sham effective potential along a particular direction, v^slab(r), would be based on the same DF as used in the calculations of the PES and describes the perturbation of the jellium. The latter would be a good approximation in particular for free-electron-like metals like Al and (somewhat less so) the noble metals Cu, Ag, Au and Ni that form part of the systems mentioned above as systems of interest for which experiments are available. If a spherically symmetrized description within a certain cut-off (muffin-tin approximation) along different directions is appropriate, phase shifts and thus scattering potential friction coefficients could be obtained analogously to eqn (3)–(5). In fact, analytical expressions have already been obtained for the anisotropic friction tensor of a diatomic molecule based on non-overlapping spherical potentials (muffin-tin approximation, eqn (C11)–(C17) of ref. 261). However, at the time (in 1975)²⁶¹ scattering phase shifts could not yet be calculated; this was first done numerically by Puska and Nieminen in 1983.¹⁶⁰ In the spirit of ref. 261 one could fit v^slab(r) to a muffin-tin form and use this to calculate the concomitant phase shifts to obtain the anisotropy and the corrugation of the friction coefficients, which can then also be compared with ODF results. Like in the LDFA, spin-polarization can be easily included in the SPF formalism, and it avoids the computational complications of ODF described above.

As noted above EF methods are not expected to be sufficient for describing the case of strong coupling, e.g. when a scattering molecule temporarily picks up an electron from the surface while strongly interacting with it. It is not yet well known under what conditions electronic friction theories break down and a strong coupling method like the IESH method should be used. This should not only depend on the value of E(CT), but also on the degree of vibrational excitation of the molecule, as stretched molecule tend to take up an electron more easily,^183,262 which is an additional reason that non-adiabatic effects are especially strong for multi-quantum vibrational relaxation.¹⁸³ Important inputs to the IESH calculations are the PESs for the diabatic neutral molecule–metal and anion–metal states, and the coupling potential. As noted, in the approach used mostly so far⁶⁵ they are obtained from an adiabatic ground state PES, the Bader charge of the molecule, and the change of the electronic energy of the system upon application of a small electric field. In an improved implementation these three quantities would all be computed with a FPB-DF fitted to first principles energies for the adiabatic ground state. As noted one can also use a pragmatic approach^169,173 in which the two diabatic PESs are computed with CDFT.^170,171 Again, an improved implementation would use a FPB-DF fitted to first principles energies for the adiabatic ground state in the CDFT calculations. Rather than using Bader charges to derive the coupling potentials as has been done so far in the CDFT approach,^169,173 CDFT can also be used to derive the coupling potential using “CDFT-couplings” (i.e., eqn (51) of ref. 171) directly, and this might be a further improvement.

3.3 PES representation and dynamics

Several approaches can be used to best represent the first principles-based results. This concerns the representation of the first principles results through the FPB density functional, as well as making an accurate fit of FPB-DFT data. In the simplest approach, one can simply fit the PES for the molecule interacting with the mobile surface using e.g. the Behler–Parinello high-dimensional neural network (BP-HDNN) method discussed above,^181,263 which uses atomic neural networks to enforce the symmetries of the system.²⁶³ Improvements can likely be made with a divide-and-conquer approach of Smits and Somers²⁶⁴ in which the full potential describing the molecule interacting with the surface is written as

V_full(r,q) = V_SS(r,q_id) + V_dist(q) + V_coup(r,q). (10)

Here, V_full(r,q) is the full interaction potential that can also be fit directly with the BP-HDNN method, and r are the molecular and q the metal atoms' coordinates. V_SS(r,q_id) is the molecule-static-surface potential with the metal atoms in their ideal lattice positions q_id. As discussed above this term can also be fitted with the corrugation reducing procedure,¹⁷⁸ or using the permutationally invariant polynomial neural network method for molecule–surface interactions.^179,180V_dist(q) is the energy required to distort the metal from its equilibrium geometry. This term can be fitted with a high-dimensional neural network method. Inspired by earlier work on QM/ME mechanical embedding,²⁶⁵ another idea²⁶⁴ one might explore is to use a highly accurate embedded atom method (EAM) fit for V_dist(q) (available for all (fcc) metal²⁶⁶ surfaces) in combination with the fits of V_coup(r,q) and of V_SS(r,q_id), with appropriate scaling of the coordinates of the metal atoms. V_coup(r,q) is the PES that describes how the interaction of the molecule with the surface changes if the lattice is distorted from its equilibrium configuration. This term can be fitted with the BP-HDNN method. By smartly combining the fits one can ensure that V_full(r,q) equals V_SS(r,q_id)for the ideal lattice configuration, and that for the molecule far away from the surface V_coup vanishes. There could be several advantages to using eqn (10). One advantage could be that an appropriately fitted V_SS(r,q_id) could be made accurate also for high interaction energies, allowing its use in TDWP calculations. Another advantage could be that the approach may well enable savings on the number of computationally expensive screened hybrid DFT calculations needed to fit a PES. Finally, we anticipate that the approach may well generate a more accurate fit of V_full(r,q) than using a brute force HDNN method to fit the full interaction directly.

The divide-and-conquer approach described above has an added advantage if it is difficult to represent accurate first principles results for differing barrier geometries with the use of a FPB density functional based on a single parameter, e.g., the parameter describing the maximum allowed EXX a. This can be mended³⁷ by allowing a to vary with impact site (X,Y) and the azimuthal orientation angle ϕ describing the molecule's orientation relative to the surface. This can be done by expanding a in a few (2–4) expansion functions of X, Y and ϕ that are totally symmetric under the plane group²⁶⁷ of the system, just like it is possible to fit a molecule-rigid surface PES using such expansion functions.²⁶⁸ An internally consistent procedure can be obtained by allowing this variation in the DFT calculations needed to fit V_SS(r,q_id) and V_coup(r,q) only. The potential describing the metal can be based on calculations using the average of a(X,Y,ϕ) over X, Y and ϕ, or the (EAM) fit already mentioned above could be used for V_dist(q). In either case a fully consistent description of the metal would be obtained. Finally, in the future it will probably be possible to derive a true first principles quality PES by adding a high-dimensional neural network (HDNN) PES based on the difference between say a thousand first principles energies and a FPB-DFT PES. This is in the spirit of the Δ-machine learning approach recently used to obtain a CCSD(T) level PES for MD simulations of liquid water,²⁶⁹ and it may already be possible to obtain DMC-quality PESs in this manner for a few selected systems of interest.

Of the systems of interest mentioned above for which experiments exist, accurate results can be obtained with the QCT method if the projectile is a diatomic molecule. For a system like D₂O + Ni(111)⁵⁹ we speculate that it might now be possible to perform TDWP calculations for the molecule interacting with the static surface, and to incorporate a posteriori corrections for the effect of surface temperature as described above in Section 2. The system incorporates nine molecular degrees of freedom. Assuming that the expense of the computation will increase by a factor of for each degree of freedom by replacing H with D, the expense of the calculation should increase relative to that performed 8 years ago for H₂O + Cu(111)²²⁰ by about a factor 2^4.5. Assuming Moore's law to hold and using that the H₂O + Cu(111)²²⁰ calculation was done 8 years ago, a full-dimensional TDWP calculation of D₂O on Ni(111) should then be possible in the year of writing (2024), ensuring that the latter system can be modeled with high accuracy. A system like CO₂ + Cu(110)²³⁶ is likely best modeled with NE-RPMD to ensure that artificial IVR between the CO stretch vibrations is avoided.¹⁹² In any case, as discussed in the previous section NE-RPMD needs more testing on systems to learn for which systems and under which conditions this method is reliable and improves significantly over QCT. A good alternative may be to use the RPH method for CO₂ + Cu(110)²³⁶ and use an a posteriori method to incorporate surface temperature. An advantage of the NE-RPMD method is that it is possible to also incorporate surface atom motion in a more direct way.¹⁹² It is also possible to incorporate electronic friction in RPMD for thermal rates,^270,271 which suggests that the same can be done for NE-RPMD, which might enable an upgraded version of MDEF with some nuclear quantum effects, like tunneling and avoidance of artificial IVR, described.

3.4 Validation and experiments

The validation of the FPB-DFT electronic structure approach in combination with the new SPF method to deal with non-adiabatic effects will require considerable care. We suggest that several systems be tackled to avoid a situation where good agreement is obtained for 1 or 2 systems for the wrong reasons, i.e., through error cancellation. The problem is that we are dealing with a situation where the ground state electronic structure method is not yet validated for systems with E(CT) < 7 eV. Specifically, DMC and RPA and the double hybrid functional XYG3 have been shown to be accurate for two systems with E(CT) > 7 eV, but we do not yet know how well they will perform for systems with E(CT) < 7 eV, and at the same time one would try to evaluate the accuracy of the new SPF method, or of another method for dealing with non-adiabatic effects. Probably the best strategy would be to test the FPB-DFT method using one specific first principles method in a systematic way on several systems, in combination with SPF, ODF, and the LDFA. If the FPB-DFT method combined with a specific EF method consistently yields the best results for all or the majority of the systems, this would constitute strong evidence for the accuracy of both the first principles electronic structure method and the particular EF method used.

The availability of experiments on both sticking and state-to-state scattering, and especially on vibrationally inelastic scattering, would be especially useful for validation of a method for describing non-adiabatic effects. The reason is as follows. Energy dissipation through, for instance, ehp excitation affects sticking and state-to-state scattering in a fundamentally different way. Whether or not a molecule undergoes DC is only affected by ehp excitation on the way to the barrier; once the molecule is over the barrier, its fate is decided (it will stick). The opposite is true for scattering: the size of e.g. state-to-state vibrationally inelastic scattering probabilities is affected by ehp excitation on the way to the barrier as well as on the way back to the gas phase. In the past this idea has been used to establish that DC of H₂ on metal surfaces is unlikely to be affected much by ehp excitation, as results for both reaction and diffractive scattering in H₂ + Pt(111) could be reproduced with an adiabatic approach using one and the same PES.²⁶² We anticipate that it could be quite useful to have accurate and detailed experiments on both sticking and vibrationally inelastic scattering available on systems like HCl + Cu(111)²³⁵ and NO + Cu(111),^237,238 as the dissociation barrier is not so high for these systems. The HCl + Au(111) system will be quite useful to work on even though the dissociation barrier is high, as experiments are available for both sticking⁴³ and scattering.^239,240 For NO on Au(111) a wealth of detailed experiments is available on scattering,^60–62,272 making this system quite useful for benchmarking, even though it might be too inert for sticking measurements. For NO + Cu(111)^237,238 some results are already available, but the accuracy of the sticking results is not so clear and not so much results are available on the vibrationally inelastic scattering yet.

We propose that barrier heights for systems with E(CT) < 7 eV, once validated and established to be accurate, are added to the systems in the SBH17 database,³² which now mostly incorporates systems with E(CT) > 7 eV (16 out of 17). The new database thus obtained would be a representative database of barriers for DC on metal surfaces. To develop a database of appropriate size, statistical methods (i.e., “least absolute shrinkage and selection operator”²⁷³ and “stepwise regression”²⁷⁴) can be employed as also used by Morgante and Peverati²⁷⁵ to reduce their large ACCDB database (8656 unique datapoints) to a much smaller database (ASCDB, 200 unique datapoints). Adding the thus obtained database to existing databases of chemical and other properties of interest,^{23,25,30,276–279} which mostly address gas phase systems, could provide considerable assistance with benchmarking density functionals that aim to be “universal”, i.e., for any system of interest to chemists and physicists. The existing large databases^{23,25,30,276–279} do contain barriers for gas phase reactions, but none of them contain barrier heights for molecule–metal surface reactions.

So far we have simply assumed that experimental data for sticking coefficients would always come from supersonic molecular beam experiments, which measure sticking probabilities at fairly well defined hyperthermal collision energies. There is also a good reason for this. Such molecular beam experiments are able to probe the reactivity on well-defined Miller index surfaces, making these experiments suitable for validation.²⁸⁰ In contrast, until recently, in the usual kinetics experiments performed under thermal conditions highly activated reactions often took place at ill-defined defects like steps, kinks, or vacancies.^281,282 This makes these older experiments unsuitable for validation purposes, because the geometry for which the barrier height needs to be computed is not clear.²⁸⁰

As nicely discussed and summarized by Jiang and co-workers,²²⁵ this situation has now changed thanks to recent experiments by Kitsopoulos, Wodtke, and Auerbach and co-workers. They have measured thermal rates for a number of reactions using a new combination of techniques, involving velocity-resolved kinetic traces determined with ion-imaging. Such experiments have addressed how steps affect CO desorption from Pt(111),²⁸³ CO desorption from Pd(111) and Pt(111),²⁸⁴ CO oxidation on Pt(111) and the stepped Pt(332) surface,^285,286 thermal CO desorption from and CO oxidation on Pd(332),²⁸⁷ NH₃ desorption from Pt(111) and Pt(332),²⁸⁸ thermal H-atom recombination on Pt(111) and Pt(332),²⁸⁹ and thermal recombination of HD on Pd(111) and Pd(332).²⁹⁰ In these experiments strategies have been employed to isolate the effects defects might have on reaction rates on macroscopic low index surfaces. For instance, rates were measured on both low index and stepped surfaces, and the assumption was made that the defect sites on the macroscopic low index surface exhibit a similar reactivity as the step sites on the stepped surface. As a result, these experiments likely offer a good testing ground to kinetics methods computing rates based on DFT modeling of the system. Inspired by these new experimental developments, new computational kinetics methods are being developed for predicting rates.^222,225,289 The CO₂ + Pd, the CO₂ + Pt and the NH₃ + Pt systems all have E(CT) < 7 eV, so the experiments on these systems are all potentially useful for testing FPB-DFT as formulated here, using appropriate computational kinetics methods for validation.

4. Summary and outlook

Dissociative chemisorption reactions on metal surfaces are of interest for both practical and scientific reasons. The barriers of these reactions are often important to the accurate modeling of heterogeneously catalyzed processes, but first principles methods capable of computing these barriers with chemical accuracy have not yet been established. Presently the molecule–metal surface interaction, which governs the molecule's reaction on as well as its state-to-state scattering from surfaces, is mostly studied with DFT. With present-day density functionals this method is not yet accurate enough, and few databases exist for testing DFT on DC barriers.

In the present state-of-the-art a semi-empirical version of DFT is used (SRP-DFT), in which a parameter in a functional with GGA exchange is fitted to reproduce sticking probabilities measured in a supersonic molecular beam experiments. Using this approach chemically accurate barriers have been extracted for a number of systems and collected in a database than can be used for testing electronic structure methods (SBH17). Unfortunately, this approach breaks down for systems affected by charge transfer, i.e., for which, as a rule of thumb, the difference between the surface's work function and the molecule's electron affinity (i.e., E(CT)) is less than 7 eV. Comparison to experiments shows that with the present approach, which is based on PESs computed with functionals containing semi-local exchange, sticking is overestimated, and vibrationally inelastic scattering is not described accurately. These systems are also prone to electronically non-adiabatic effects like electron–hole pair excitations, for which the accuracy of existing methods is difficult if not impossible to benchmark separately.

Because there are two sources of uncertainty in computing sticking probabilities for systems with E(CT) < 7 eV, a semi-empirical approach for adjusting the functional to only reproduce measured sticking probabilities cannot be expected to work. This is unfortunate because these systems are quite important to sustainable chemistry. To address the problem described above, this Perspective has first summarized the present state-of-the-art in electronic structure theory for molecules interacting with metal surfaces, in the description of non-adiabatic effects in these systems, and in extracting computed observables for comparison with experiments through fitting of PESs and dynamics calculations. After that a vision has been sketched on how to make progress in these areas so as to eventually achieve accurate theoretical predictions for systems that are often of high practical relevance to sustainable chemistry.

Concerning electronic structure theory, a problem with DFT using standard semi-local functionals (GGAs or meta-GGAs) is that such functionals are usually inaccurate for predicting sticking curves. This implies that they are inaccurate for barrier heights to DC on metals. For systems with E(CT) > 7 eV this can be addressed by constructing system-specific parameterized functionals, which mix semi-local exchange functionals and use either semi-local or non-local but efficient-to-evaluate correlation functionals. In this semi-empirical approach, the mixing parameter is varied until sticking probabilities measured in supersonic molecular beam experiments are reproduced. For these systems this approach works because the minimum barrier height can be straddled with functionals containing semi-local exchange, and because standard functionals are capable of describing how the barrier height to sticking varies with system geometry.

The semi-empirical approach described above does not work for systems with E(CT) < 7 eV because functionals using semi-local exchange systematically underestimate barrier heights for these systems. Calculations on an infamous system with E(CT) < 7 eV, i.e., O₂ + Al(111), suggest that this problem may be resolved by using screened hybrid density functionals. We have also noted a fundamental DFT problem with describing molecules interacting with metals: at long range the fraction of exact exchange should be maximum in the gas phase, while it should be minimum (exact exchange should be screened) in the metal.

Two first principles, or non-empirical, methods have recently been demonstrated to show promise of predictive accuracy for barrier heights in two benchmark systems. Diffusion Monte-Carlo achieved an accuracy of about 1.5 kcal mol⁻¹ for H₂ + Cu(111) and Al(110). The RPA method showed chemical accuracy for both systems (errors <1 kcal mol⁻¹). Because results are only available for two systems at this stage, it cannot yet be said which method will be best; DMC worked better than the RPA for a database with 76 barrier heights for gas phase reactions.

An application of the CCSD(T) method, which has been called the gold standard for computing gas phase barriers, to DC on metals has not yet been demonstrated. Like other ab initio many electron-wave function methods (or correlated wave function methods), CCSD(T) is currently computationally too expensive for periodic calculations on DC on metals. However, such methods show promise when applied in an embedded cluster fashion, using density functional embedding. But the accuracy of this embedded correlated wave function approach is not yet as high as achieved with the RPA and DMC methods. Another cluster-based approach, the ONIOM method, has shown very promising results when carefully monitoring the convergence with respect to cluster size and combined with a doubly hybrid functional. Like the RPA, the latter makes the PES depend on unoccupied states and thus comes at a computational price that appears to be unavoidable when high accuracy is needed.

Systems in which a molecule approaches a metal surface are always prone to electron–hole pair excitation due to molecular motion, an electronically non-adiabatic effect breaking the Born–Oppenheimer approximation. If the coupling between molecular motion and the electrons in the metal is weak, this can be addressed with electronic friction theory. Friction tensors can be incorporated in a generalized Langevin form for the nuclear equations of motion. Two EF methods currently exist. The ODF method takes the electronic structure of the molecule and the metal surface into account. Problems of the ODF method are that it exhibits a non-physical dependence on the broadening parameters needed to compute friction coefficients from electron–phonon coupling matrix elements, and that its use may lead to unphysically large friction in regions where a spin transition of an impinging atom or molecule occurs. Going beyond the Markov approximation that is currently applied to coarse grain the effect of the electron on the nuclear dynamics could be a way forward. The LDFA method does not exhibit these problems, and has been shown to be reliable for atoms scattering from metals. However, the LDFA method does not include effects related to the electronic structure of the molecule and of the metal other than the mere perturbation of an atom embedded in (bulk) jellium. Neither of these two methods has at present been proven to be universally more accurate for describing non-adiabatic effects on DC on metal surfaces. In fact, they give quite different results for at least one benchmark system, i.e., N₂ + Ru(0001).

Cases in which the coupling between molecular motion and the metal electrons is strong, for instance when the molecule is able to (temporarily) pick up a partial charge from the metal in what is not the system's electronic ground state, cannot be described by electronic friction methods. Such cases can be described with the IESH method, in which the motion of the neutral molecule on the metal surface is coupled to states in which the molecule is an anion or the metal electron is excited to a virtual metal level. In the present state-of-the-art the two molecule–metal surface states are computed with CDFT, and the coupling potentials with a model involving Bader charges and the change of energy of the charged system in an imposed electric field.

Concerning dynamics to compute observables like sticking coefficients and scattering probabilities, the fitting of a PES can be avoided in direct dynamics calculations. This is feasible only if the electronic structure approach is computationally efficient and if the probabilities to be computed are large enough that they can be computed with high enough statistical accuracy using a limited number of classical trajectories. Otherwise a divide-and-conquer approach needs to be used in which electronic structure data are computed first and next fitted to a PES, after which calculations using an appropriate dynamical model and method need to be performed. For PES fitting accurate methods are now available, and this stage does usually not present bottlenecks if enough electronic structure data are available for making an accurate fit.

In the dynamical model, a rule of thumb is that surface atom motion has to be modeled for molecules heavier than D₂, and/or for surface temperatures considerably higher than room temperature. It is advisable to attempt to model non-adiabatic effects for systems for with E(CT) < 7 eV. Both (nuclear and electronic) dissipative degrees of freedom can be modeled computing forces on the fly with AIMDEF or with a pre-computed PES using MDEF.

Of the dynamics methods, the QCT method is usually highly accurate at describing activated sticking in supersonic molecular beam experiments. The QCT method may also yield accurate results for sticking if improved methods are used for assigning final states of scattered molecules and the so-called adiabatic correction is applied. In specific cases (low nozzle temperature, or the vibrationally excited state involves the highest frequency vibration of the molecule isolated from other vibrations) the QCT method may also yield accurate results for sticking of polyatomic molecules.

Concerning quantum dynamics, the TDWP method is a very accurate method that has been applied to DC of diatomic molecules and of the H₂O molecule. The computational expense of the method is currently too high to treat sticking of bigger molecules. Sticking of intermediately sized molecules can be modeled quantum dynamically with the reaction path Hamiltonian method. With both the TDWP and the RPH method it is possible to rather accurately describe the effects of surface temperature using post-processing methods. Researchers are also starting to explore the accuracy of non-equilibrium ring polymer molecular dynamics, which can describe sticking in the tunneling regime more accurately than QCT. The use of NE-RPMD avoids artificial intra-molecular vibrational relaxation in the molecule on the way to the surface, which is a known problem with the QCT method. However, NE-RPMD is not more accurate than QCT under all conditions, and more research is needed to establish empirically under which conditions and for which systems NE-RPMD improves over the QCT method.

Coming to the way forward, it is clear that the greatest challenges exist for systems with E(CT) < 7 eV. A number of systems have been identified for which accurate results are available from supersonic molecular beam experiments, or for which accurate rates are available from thermal experiments using new techniques (velocity-resolved kinetic traces determined with ion imaging). Calculations with the new approaches suggested in this Perspective can model these systems for validation.

The greatest challenge likely lies with the electronic structure approach for the electronic ground state. As a steppingstone for testing first principles-based methods, we suggest to use system-specific parameterized density functionals as before, but now to base the parameters in these functionals on calculations with accurate first principles methods, like DMC or the RPA method (FPB-DFT), for judiciously chosen geometries. Systematic ways of diminishing the fixed-node and the locality errors in DMC exist. It will probably be best to use a parameterized density functional with exchange taken from a screened hybrid functional, combined with one of the Chalmers–Rutgers van de Waals correlation functionals. Judicious choices exist for the combination of these, and the maximum fraction of exchange and/or the range parameter in the functional can be tuned to fit to the first principles result.

The next great challenge is to come up with an EF method that combines the advantages and avoids the pitfalls of the two existing methods, the LDFA and the ODF methods. With the new (SPF) approach, we suggest to focus on extracting an electronic scattering potential (as a Kohn–Sham effective potential) from a DFT calculation for the full molecule–metal surface system. The method would use analytical expressions that were already derived, but at a time when it was not yet possible to evaluate them numerically on computers. Improvements of the IESH approach have also been discussed and consist of using a FPB density functional in the CDFT calculations, and actual CDFT couplings for the non-adiabatic coupling potentials.

For fitting PESs we suggest using the high-dimensional neural network potential approach of Behler and Parinello, or more recently developed machine learning techniques with different descriptor representations. Improvements can likely be made by splitting the molecule–surface interaction up in three components. In this way, if needed one might also correct for errors in fitting the density functional to first principles results at different geometries, by allowing the fitting parameter to vary according to the full symmetry of the molecule–metal surface system. It may also be possible to upgrade PESs obtained with FPB-DFT to an actual first principles PES. This can likely be done using a Δ-machine learning approach in which, for a limited number of judiciously chosen points, an added neural network potential is fit to the difference between first principles and fitted FPB-DFT results.

We have also provided some discussion on how some of the experimental results we suggest to model can be best addressed with available dynamics methods. For instance, we suggest that the D₂O + Ni(111) system can now be tackled with the TDWP method, and that calculations on CO₂ + Cu(110) can be done with NE-RPMD, and/or the RPH method.

Furthermore, we have discussed ways of validating the new computational approaches with comparisons to experiments. We have argued that systems for which experimental sticking and vibrational state-to-state scattering results exist will be especially useful to validating approaches for electronically non-adiabatic scattering. The reason is that these two processes are affected by non-adiabatic energy dissipation in fundamentally different ways.

Once accurate barrier heights are available for systems with E(CT) < 7 eV, these should be added to results that are already available for systems with E(CT) > 7 eV in the SBH17 database. This way, the first representative (i.e., including systems with no restrictions on E(CT)) database for barrier heights to DC on metal surfaces would be obtained. This database would be quite useful to computational heterogeneous catalysis. Adding this database to existing large databases containing results for mostly gas phase systems should also be very useful, as it should clearly be desirable to be able to test whether new electronic structure approaches work for both types of systems. Such a database would also be useful to testing a new “true made simple hybrid density functional”, which would incorporate the correct limiting behavior of long-range exact exchange in both the gas phase and the metal. We anticipate that such a functional should be a local hybrid, using diagnostics based on the kinetic energy density to determine where the electrons are in the system to more accurately evaluate the exchange interaction between them. We also anticipate that with the right training machine-learning based density functionals incorporating exchange energy densities might be successful at exhibiting the right long-range behavior of exchange in both parts of molecule–metal surface systems.

In summary, a clear route can be envisaged to an improved modeling of molecule–metal surface systems that are prone to charge transfer, many of which are likely important to sustainable chemistry. As always, the devil will be in the details and, as always, it will be fun coming up with the solutions with well thought out fundamental research. We anticipate this to be a vibrant direction of research that can provide wonderful challenges to researchers for decades to come before they can declare “problem solved” and go on their way to solve the next great problem.

Author contributions

G. J. K. wrote the outline of the manuscript. Both G. J. K. and J. M. contributed to the writing of the original draft of the manuscript, the review of the manuscript, and the editing of the manuscript.

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

We are grateful to Prof. Bin Jiang for providing us with an electronic copy of Fig. 12. Part of the work presented here was supported financially by the ERC (with the ERC-ADG-2013 grant no. 338580) and by NWO (through CW-TOP grants no. 715.011.000 and 715.017.001) and through several grants of computer time by NWO-NCF and NWO-EW.

References

1. R. Noyori, Nat. Chem., 2009, 1, 5–6.
2. C. A. Wolcott, A. J. Medford, F. Studt and C. T. Campbell, J. Catal., 2015, 330, 197–207.
3. M. K. Sabbe, M. F. Reyniers and K. Reuter, Catal. Sci. Technol., 2012, 2, 2010–2024.
4. R. D. Beck, P. Maroni, D. C. Papageorgopoulos, T. T. Dang, M. P. Schmid and T. R. Rizzo, Science, 2003, 302, 98–100.
5. G. Ertl, J. Vac. Sci. Technol., A, 1983, 1, 1247–1253.
6. K. Honkala, A. Hellman, I. N. Remediakis, A. Logadottir, A. Carlsson, S. Dahl, C. H. Cristensen and J. K. Nørskov, Science, 2005, 307, 555–558.
7. A. Bruix, J. T. Margraf, M. Andersen and K. Reuter, Nat. Catal., 2019, 2, 659–670.
8. B. W. J. Chen, L. Xu and M. Mavrikakis, Chem. Rev., 2021, 121, 1007–1048.
9. A. J. Medford, A. Vojvodic, J. S. Hummelshøj, J. Voss, F. Abild-Pedersen, F. Studt, T. Bligaard, A. Nilsson and J. K. Nørskov, J. Catal., 2015, 328, 36–42.
10. A. J. Medford, J. Wellendorff, A. Vojvodic, F. Studt, F. Abild-Pedersen, K. W. Jacobsen, T. Bligaard and J. K. Nørskov, Science, 2014, 345, 197–200.
11. J. K. Nørskov, F. Abild-Petersen, F. Studt and T. Bligaard, Proc. Natl. Acad. Sci. U.S.A., 2011, 108, 937–943.
12. S. Perathoner, G. Centi, S. Gross and E. M. J. Hensen, Science and technology roadmap on catalysis for Europa: A path to create a sustainable future, https://www.euchems.eu/wp-content/uploads/2016/07/160729-Science-and-Technology-Roadmap-on-Catalysis-for-Europe-2016.pdf.
13. G. J. Kroes, Phys. Chem. Chem. Phys., 2021, 23, 8962–9048.
14. B. Jiang and H. Guo, J. Chem. Phys., 2019, 150, 180901.
15. K. Golibrzuch, N. Bartels, D. J. Auerbach and A. M. Wodtke, Annu. Rev. Phys. Chem., 2015, 66, 399–425.
16. D. Farías and R. Miranda, Prog. Surf. Sci., 2011, 86, 222–254.
17. H. Chadwick and R. D. Beck, Annu. Rev. Phys. Chem., 2017, 68, 39–61.
18. D. J. Auerbach, J. C. Tully and A. M. Wodtke, Nat. Sci., 2021, 1, e10005.
19. M. Alducin, R. Diéz Muiño and J. I. Juaristi, Prog. Surf. Sci., 2017, 92, 317–340.
20. J. P. Perdew, MRS Bull., 2013, 38, 743–750.
21. A. J. Cohen, P. Mori-Sánchez and W. T. Yang, Science, 2008, 321, 792–794.
22. E. Sim, S. Song, S. Vuckovic and K. Burke, J. Am. Chem. Soc., 2022, 144, 6625–6639.
23. N. Mardirossian and M. Head-Gordon, Mol. Phys., 2017, 115, 2315–2372.
24. R. Peverati and D. G. Truhlar, Philos. Trans. R. Soc., A, 2014, 372, 20120476.
25. L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 32184–32214.
26. F. Libisch, C. Huang, P. L. Liao, M. Pavone and E. A. Carter, Phys. Rev. Lett., 2012, 109, 198303.
27. K. Doblhoff-Dier, J. Meyer, P. E. Hoggan and G. J. Kroes, J. Chem. Theory Comput., 2017, 13, 3208–3219.
28. Z. Y. Wei, J. M. P. Martirez and E. A. Carter, J. Chem. Phys., 2023, 159, 194108.
29. A. Pribram-Jones, D. A. Gross and K. Burke, Annu. Rev. Phys. Chem., 2015, 66, 283–304.
30. P. Morgante and R. Peverati, J. Comput. Chem., 2019, 40, 839–848.
31. S. Mallikarjun Sharada, T. Bligaard, A. C. Luntz, G. J. Kroes and J. K. Nørskov, J. Phys. Chem. C, 2017, 121, 19807–19815.
32. T. Tchakoua, N. Gerrits, E. W. F. Smeets and G. J. Kroes, J. Chem. Theory Comput., 2023, 19, 245–270.
33. C. Díaz, E. Pijper, R. A. Olsen, H. F. Busnengo, D. J. Auerbach and G. J. Kroes, Science, 2009, 326, 832–834.
34. F. Nattino, D. Migliorini, G. J. Kroes, E. Dombrowski, E. A. High, D. R. Killelea and A. L. Utz, J. Phys. Chem. Lett., 2016, 7, 2402–2406.
35. D. Migliorini, H. Chadwick, F. Nattino, A. Gutiérrez-Gonzalez, E. Dombrowski, E. A. High, H. Guo, A. L. Utz, B. Jackson, R. D. Beck and G. J. Kroes, J. Phys. Chem. Lett., 2017, 8, 4177–4182.
36. T. Tchakoua, T. Jansen, Y. van Nies, R. F. A. van den Elshout, B. A. B. van Boxmeer, S. P. Poort, M. G. Ackermans, G. S. Beltrão, S. A. Hildebrand, S. E. J. Beekman, T. van der Drift, S. Kaart, A. Santic, E. E. Spuijbroeck, N. Gerrits, M. F. Somers and G. J. Kroes, J. Phys. Chem. A, 2023, 127, 10481–10498.
37. A. D. Powell, N. Gerrits, T. Tchakoua, M. F. Somers, H. F. Busnengo, J. Meyer, G. J. Kroes and K. Doblhoff-Dier, J. Phys. Chem. Lett., 2024, 15, 307–315.
38. A. D. Powell, G. J. Kroes and K. Doblhoff-Dier, J. Chem. Phys., 2020, 153, 224701.
39. M. Karikorpi, S. Holloway, N. Henriksen and J. K. Nørskov, Surf. Sci., 1987, 179, L41–L48.
40. A. Salin, J. Chem. Phys., 2006, 124, 104704.
41. N. Gerrits, E. W. F. Smeets, S. Vuckovic, A. D. Powell, K. Doblhoff-Dier and G. J. Kroes, J. Phys. Chem. Lett., 2020, 11, 10552–10560.
42. K. R. Bryenton, A. A. Adeleke, S. G. Dale and E. R. Johnson, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2022, 13, e1631.
43. N. Gerrits, J. Geweke, E. W. F. Smeets, J. Voss, A. M. Wodtke and G. J. Kroes, J. Phys. Chem. C, 2020, 124, 15944–15960.
44. Q. H. Liu, X. Y. Zhou, L. S. Zhou, Y. L. Zhang, X. Luo, J. Guo and B. Jiang, J. Phys. Chem. C, 2018, 122, 1761–1769.
45. G. Füchsel, X. Y. Zhou, B. Jiang, J. I. Juaristi, M. Alducin, H. Guo and G. J. Kroes, J. Phys. Chem. C, 2019, 123, 2287–2299.
46. P. R. Shirhatti, J. Geweke, C. Steinsiek, C. Bartels, I. Rahinov, D. J. Auerbach and A. M. Wodtke, J. Phys. Chem. Lett., 2016, 7, 1346–1350.
47. T. H. Liu, B. N. Fu and D. H. Zhang, J. Chem. Phys., 2017, 140, 164706.
48. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413–7421.
49. Q. Shen, L. Zhu, J. Wu, W. Dong, X. Wang, T. Wang, B. Jiang and X. M. Yang, Chin. J. Chem. Phys., 2024, 37, 490–496.
50. N. Gerrits, E. W. F. Smeets, S. Vuckovic, A. D. Powell, K. Doblhoff-Dier and G. J. Kroes, J. Phys. Chem. Lett., 2020, 11, 10552–10560.
51. J. Behler, B. Delley, S. Lorenz, K. Reuter and M. Scheffler, Phys. Rev. Lett., 2005, 94, 036104.
52. J. Behler, K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 115421.
53. M. Ramos, C. Díaz, A. E. Martínez, H. F. Busnengo and F. Martín, Phys. Chem. Chem. Phys., 2017, 19, 10217–10221.
54. B. Jiang and H. Guo, Phys. Chem. Chem. Phys., 2016, 18, 21817–21824.
55. X. X. Hu, M. Yang, D. Q. Xi and H. Guo, J. Chem. Phys., 2018, 149, 044703.
56. N. Gerrits and G. J. Kroes, J. Phys. Chem. C, 2019, 123, 28291–28300.
57. M. Pozzo and D. Alfè, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 245313.
58. B. Oudot and K. Doblhoff-Dier, J. Chem. Phys., 2024, 161, 054708.
59. P. M. Hundt, B. Jiang, M. E. van Reijzen, H. Guo and R. D. Beck, Science, 2014, 344, 504–507.
60. Y. Huang, A. M. Wodtke, H. Hou, C. T. Rettner and D. J. Auerbach, Phys. Rev. Lett., 2000, 84, 2985–2988.
61. K. Golibrzuch, P. R. Shirhatti, I. Rahinov, A. Kandratsenka, D. J. Auerbach, A. M. Wodtke and C. Bartels, J. Chem. Phys., 2014, 140, 044701.
62. N. Bartels, B. C. Krüger, D. J. Auerbach, A. M. Wodtke and T. Schäfer, Angew. Chem., Int. Ed., 2014, 53, 13690–13694.
63. C. L. Box, Y. L. Zhang, R. R. Yin, B. Jiang and R. J. Maurer, JACS Au, 2021, 1, 164–173.
64. R. R. Yin, Y. L. Zhang and B. Jiang, J. Phys. Chem. Lett., 2019, 10, 5969–5974.
65. S. Roy, N. Shenvi and J. C. Tully, J. Chem. Phys., 2009, 130, 174716.
66. M. Head-Gordon and J. C. Tully, J. Chem. Phys., 1995, 103, 10137–10145.
67. N. Shenvi, S. Roy and J. C. Tully, Science, 2009, 326, 829–832.
68. N. Shenvi, S. Roy and J. C. Tully, J. Chem. Phys., 2009, 130, 174107.
69. P. Spiering, K. Shakouri, J. Behler, G. J. Kroes and J. Meyer, J. Phys. Chem. Lett., 2019, 10, 2957–2962.
70. K. Shakouri, J. Behler, J. Meyer and G. J. Kroes, J. Phys. Chem. Lett., 2017, 8, 2131–2136.
71. M. M. Montemore, M. A. van Spronsen, R. J. Madix and C. M. Friend, Chem. Rev., 2017, 118, 2816–2862.
72. A. Morales-Garcia, F. Viñes, J. R. B. Gomes and F. Illas, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2021, 11, e1530.
73. D. H. Zhang and H. Guo, Annu. Rev. Phys. Chem., 2016, 67, 135–158.
74. T. M. Henderson, A. F. Izmaylov, G. E. Scuseria and A. Savin, J. Chem. Theory Comput., 2008, 4, 1254–1262.
75. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
76. K. Lee, É. D. Murray, L. Z. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081101.
77. R. Sabatini, T. Gorni and S. de Gironcoli, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 041108.
78. H. A. Michelsen, C. T. Rettner, D. J. Auerbach and R. N. Zare, J. Chem. Phys., 1993, 98, 8294–8307.
79. E. W. F. Smeets, J. Voss and G. J. Kroes, J. Phys. Chem. A, 2019, 123, 5395–5406.
80. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
81. J. W. Sun, A. Ruzsinszky and J. P. Perdew, Phys. Rev. Lett., 2015, 115, 036402.
82. J. M. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
83. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Constantin and J. W. Sun, Phys. Rev. Lett., 2009, 103, 026403.
84. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. L. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
85. E. W. F. Smeets and G. J. Kroes, J. Phys. Chem. C, 2021, 125, 8993–9010.
86. E. N. Ghassemi, M. Wijzenbroek, M. F. Somers and G. J. Kroes, Chem. Phys. Lett., 2017, 683, 329–335.
87. A. C. Luntz, J. K. Brown and M. D. Williams, J. Chem. Phys., 1990, 93, 5240–5246.
88. J. W. Sun, R. Haunschild, B. Xiao, I. W. Bulik, G. E. Scuseria and J. P. Perdew, J. Chem. Phys., 2013, 138, 044113.
89. J. Wellendorff, K. T. Lundgaard, A. Møgelhøj, V. Petzold, D. D. Landis, J. K. Nørskov, T. Bligaard and K. W. Jacobsen, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 235149.
90. B. J. Lynch, P. L. Fast, M. Harris and D. G. Truhlar, J. Phys. Chem. A, 2000, 104, 4811–4815.
91. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
92. W. W. Gao, T. A. Abtew, T. Cai, Y. Y. Sun, S. B. Zhang and P. H. Zhang, Solid State Commun., 2016, 234–235, 10–13.
93. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
94. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2006, 124, 219906.
95. R. R. Yin, Y. L. Zhang, F. Libisch, E. A. Carter, H. Guo and B. Jiang, J. Phys. Chem. Lett., 2018, 9, 3271–3277.
96. L. Österlund, I. Zoric and B. Kasemo, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 15452–15455.
97. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
98. R. A. B. van Bree, N. Gerrits and G. J. Kroes, Faraday Discuss., 2024, 251, 361–381.
99. D. Bohm and D. Pines, Phys. Rev., 1951, 82, 625–634.
100. D. Pines and D. Bohm, Phys. Rev., 1952, 85, 338–353.
101. D. Bohm and D. Pines, Phys. Rev., 1953, 92, 609–625.
102. D. C. Langreth and J. P. Perdew, Phys. Rev. B: Solid State, 1977, 15, 2884–2901.
103. J. Paier, X. Ren, P. Rinke, G. E. Scuseria, A. Grüneis, G. Kresse and M. Scheffler, New J. Phys., 2012, 14, 043002.
104. X. J. Zhou and F. Wang, J. Comput. Chem., 2017, 38, 798–806.
105. K. Krongchon, B. Busemeyer and L. K. Wagner, J. Chem. Phys., 2017, 146, 124129.
106. J. Harl, L. Schimka and G. Kresse, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 115126.
107. Y. Zhang, X. Xu and W. A. G. III, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 4963–4968.
108. Z. Chen, Z. Liu and X. Xu, Nat. Commun., 2023, 14, 936.
109. N. Mardirossian and M. Head-Gordon, J. Chem. Phys., 2016, 144, 214110.
110. C.-O. Almbladh and U. von Barth, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 3231–3244.
111. T. M. Henderson, A. F. Izmaylov, G. E. Scuseria and A. Savin, J. Chem. Phys., 2007, 127, 221103.
112. B. M. Austin, D. Y. Zubarev and W. A. Lester, Jr., Chem. Rev., 2012, 112, 263–288.
113. P. E. Hoggan, Adv. Quantum Chem., 2018, 76, 271–278.
114. R. O. Sharma, T. T. Rantala and P. E. Hoggan, J. Phys. Chem. C, 2020, 124, 26232–26240.
115. K. Raghavachari, G. W. Trucks, J. A. Pople and M. Headgordon, Chem. Phys. Lett., 1989, 157, 479–483.
116. A. Kubas, D. Berger, H. Oberhofer, D. Maganas, K. Reuter and F. Neese, J. Phys. Chem. Lett., 2016, 7, 4207–4212.
117. H. Z. Ye and T. C. Berkelbach, Faraday Discuss., 2024, 254, 628–640.
118. B. X. Shi, A. Zen, V. Kapil, P. R. Nagy, A. Grüneis and A. Michaelides, J. Am. Chem. Soc., 2023, 145, 25372–25381.
119. T. Tsatsoulis, S. Sakong, A. Groß and A. Grüneis, J. Chem. Phys., 2018, 149, 244105.
120. I. Y. Zhang and A. Grüneis, Front. Mater., 2019, 6, 123.
121. F. Libisch, C. Huang and E. A. Carter, Acc. Chem. Res., 2014, 47, 2768–2775.
122. K. Andersson, P. Å. Malmqvist, B. O. Roos, A. J. Sadlej and K. Wolinski, J. Phys. Chem., 1990, 94, 5483–5488.
123. P. Celani and H. J. Werner, J. Chem. Phys., 2000, 112, 5546–5557.
124. Q. Zhao, X. Zhang, J. M. P. Martirez and E. A. Carter, J. Chem. Theory Comput., 2020, 16, 7078–7088.
125. C. T. Rettner, H. A. Michelsen and D. J. Auerbach, J. Chem. Phys., 1995, 102, 4625–4641.
126. C. Angeli, R. Cirimaglia, S. Evangelisti, T. Leininger and J. P. Malrieu, J. Chem. Phys., 2001, 114, 10252–10264.
127. K. Morokuma, Bull. Korean Chem. Soc., 2003, 24, 797–801.
128. F. Göltl, C. Houriez, M. Guitou, G. Chambaud and P. Sautet, J. Phys. Chem. C, 2014, 118, 5374–5382.
129. H. J. Werner and P. J. Knowles, J. Chem. Phys., 1988, 89, 5803–5814.
130. S. R. Langhoff and E. R. Davidson, Int. J. Quantum Chem., 1974, 8, 61–72.
131. R. B. Araujo, G. L. S. Rodrigues, E. Campos dos Santos and L. G. M. Petterson, Nat. Commun., 2022, 13, 6853.
132. S. Grimme, J. Anthony, S. Ehrlich and K. Goerigk, J. Chem. Phys., 2010, 132, 154104.
133. A. Allouche, J. Comput. Chem., 2012, 32, 174–182.
134. E. Hasselbrink, Science, 2009, 326, 809–810.
135. E. Müller-Hartmann, T. V. Ramakrishnan and G. Toulouse, Phys. Rev. B: Solid State, 1971, 3, 1102–1119.
136. S. P. Rittmeyer, V. J. Bukas and K. Reuter, Adv. Phys.: X, 2018, 3, 1381574.
137. I. G. Ryabinkin and A. F. Izmaylov, J. Phys. Chem. Lett., 2017, 8, 440–444.
138. Z. Jin and J. E. Subotnik, J. Chem. Theory Comput., 2021, 17, 614–626.
139. J. C. Tully, J. Chem. Phys., 2024, 160, 124708.
140. Y. Li and G. Wahnström, Phys. Rev. Lett., 1992, 68, 3444–3447.
141. W. Dou and J. E. Subotnik, J. Chem. Phys., 2018, 148, 230901.
142. S. P. Rittmeyer, J. Meyer and K. Reuter, Phys. Rev. Lett., 2017, 119, 176808.
143. B. Hellsing and M. Persson, Phys. Scr., 1984, 29, 360–371.
144. M. Askerka, R. J. Maurer, V. S. Batista and J. C. Tully, Phys. Rev. Lett., 2016, 116, 217601.
145. R. J. Maurer, M. Askerka, V. S. Batista and J. C. Tully, Phys. Rev. B, 2016, 94, 115432.
146. P. Spiering and J. Meyer, J. Phys. Chem. Lett., 2018, 9, 1803–1808.
147. C. L. Box, W. G. Stark and R. J. Maurer, Electron. Struct., 2023, 5, 035005.
148. Y. L. Zhang, R. J. Maurer and B. Jiang, J. Phys. Chem. C, 2020, 124, 186–195.
149. M. Persson and B. Hellsing, Phys. Rev. Lett., 1982, 49, 662–665.
150. J. R. Trail, M. C. Graham and D. M. Bird, Comput. Phys. Commun., 2001, 137, 163–173.
151. D. Novko, M. Alducin, M. Blanco-Rey and J. I. Juaristi, Phys. Rev. B, 2016, 94, 224306.
152. Y. L. Zhang, R. J. Maurer, H. Guo and B. Jiang, Chem. Sci., 2019, 10, 1089–1097.
153. R. J. Maurer, B. Jiang, H. Guo and J. C. Tully, Phys. Rev. Lett., 2017, 118, 256001.
154. J. R. Trail, D. M. Bird, M. Persson and S. Holloway, J. Chem. Phys., 2003, 119, 4539–4549.
155. J. R. Trail, M. C. Graham, D. M. Bird, M. Persson and S. Holloway, Phys. Rev. Lett., 2002, 88, 166802.
156. A. C. Luntz and M. Persson, J. Chem. Phys., 2005, 123, 074704.
157. A. C. Luntz, M. Persson and G. O. Sitz, J. Chem. Phys., 2006, 124, 091101.
158. G. Füchsel, T. Klamroth, S. Monturet and P. Saalfrank, Phys. Chem. Chem. Phys., 2011, 13, 8659–8670.
159. J. I. Juaristi, M. Alducin, R. Díez Muiño, H. F. Busnengo and A. Salin, Phys. Rev. Lett., 2008, 100, 116102.
160. M. J. Puska and R. M. Nieminen, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 27, 6121–6128.
161. N. Gerrits, J. I. Juaristi and J. Meyer, Phys. Rev. B, 2020, 102, 155130.
162. A. C. Luntz, I. Makkonen, M. Persson, S. Holloway, D. M. Bird and M. S. Mizielinski, Phys. Rev. Lett., 2009, 102, 109601.
163. P. M. Echenique, R. M. Nieminen and R. H. Ritchie, Solid State Commun., 1981, 37, 779–781.
164. A. Salin, A. Arnau, P. M. Echenique and E. Zaremba, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 2537–2548.
165. E. W. F. Smeets, G. Füchsel and G. J. Kroes, J. Phys. Chem. C, 2019, 123, 23049–23063.
166. E. N. Ghassemi, E. W. F. Smeets, M. F. Somers, G. J. Kroes, I. M. N. Groot, L. B. F. Juurlink and G. Füchsel, J. Phys. Chem. C, 2019, 123, 2973–2986.
167. S. Roy, N. Shenvi and J. C. Tully, J. Phys. Chem. C, 2009, 113, 16311–16320.
168. J. C. Tully, J. Chem. Phys., 1990, 93, 1061–1071.
169. G. Meng and B. Jiang, J. Chem. Phys., 2022, 157, 214103.
170. Q. Wu and T. Van Voorhis, J. Chem. Theory Comput., 2006, 2, 765–774.
171. B. Kaduk, T. Kowalczyk and T. Van Voorhis, Chem. Rev., 2011, 112, 321–370.
172. Y. Zhang and W. Yang, Phys. Rev. Lett., 1998, 80, 890.
173. G. Meng, J. Gardner, N. Hertl, W. Dou, R. J. Maurer and B. Jiang, Phys. Rev. Lett., 2024, 133, 036203.
174. R. J. V. Wagner, N. Henning, B. C. Krüger, G. B. Park, J. Altschaffel, A. Kandratsenka, A. M. Wodtke and T. Schafer, J. Phys. Chem. Lett., 2017, 8, 4887–4892.
175. C. Díaz, R. A. Olsen, H. F. Busnengo and G. J. Kroes, J. Phys. Chem. C, 2010, 114, 11192–11201.
176. A. Groß and A. Dianat, Phys. Rev. Lett., 2007, 98, 206107.
177. X. Y. Zhou, B. Kolb, X. Luo, H. Guo and B. Jiang, J. Phys. Chem. C, 2017, 121, 5594–5602.
178. H. F. Busnengo, A. Salin and W. Dong, J. Chem. Phys., 2000, 112, 7641–7651.
179. B. Jiang and H. Guo, J. Chem. Phys., 2014, 141, 034109.
180. B. Jiang, J. Li and H. Guo, Int. Rev. Phys. Chem., 2016, 35, 479–506.
181. J. Behler and M. Parrinello, Phys. Rev. Lett., 2007, 98, 146401.
182. Y. L. Zhang, C. Hu and B. Jiang, J. Phys. Chem. Lett., 2019, 10, 4962–4967.
183. Y. H. Huang, C. T. Rettner, D. J. Auerbach and A. M. Wodtke, Science, 2000, 290, 111–114.
184. A. M. Wodtke, D. Matsiev and D. J. Auerbach, Prog. Surf. Sci., 2008, 83, 167–214.
185. K. Shakouri, J. Behler, J. Meyer and G. J. Kroes, J. Phys. Chem. C, 2018, 122, 23470–23480.
186. M. Blanco-Rey, J. I. Juaristi, R. Díez Muiño, H. F. Busnengo, G. J. Kroes and M. Alducin, Phys. Rev. Lett., 2014, 112, 103203.
187. P. Saalfrank, J. I. Juaristi, M. Alducin, M. Blanco-Rey and R. Diéz Muiño, J. Chem. Phys., 2014, 141, 234702.
188. R. N. Porter and L. M. Raff, in Dynamics of Molecular Collisions, Part B, ed. W. H. Miller, Plenum, New York, 1976, pp. 1–52.
189. E. N. Ghassemi, M. Somers and G. J. Kroes, J. Phys. Chem. C, 2018, 122, 22939–22952.
190. T. Tchakoua, A. D. Powell, N. Gerrits, M. F. Somers, K. Doblhoff-Dier, H. F. Busnengo and G. J. Kroes, J. Phys. Chem. C, 2023, 127, 5395–5407.
191. F. Nattino, H. Ueta, H. Chadwick, M. E. van Reijzen, R. D. Beck, B. Jackson, M. C. van Hemert and G. J. Kroes, J. Phys. Chem. Lett., 2014, 5, 1294–1299.
192. N. Gerrits, B. Jackson and A. Bogaerts, J. Phys. Chem. Lett., 2024, 15, 2566–2572.
193. F. Nattino, D. Migliorini, M. Bonfanti and G. J. Kroes, J. Chem. Phys., 2016, 144, 044702.
194. L. Bonnet and J. C. Rayez, Chem. Phys. Lett., 1997, 277, 183–190.
195. L. Bonnet, Int. Rev. Phys. Chem., 2013, 32, 171–228.
196. A. Rodríguez-Fernández, L. Bonnet, C. Crespos, P. Larrégaray and R. Diéz Muiño, J. Phys. Chem. Lett., 2019, 10, 7629–7635.
197. L. Zhang, J. Chen and B. Jiang, J. Phys. Chem. C, 2021, 125, 4995–5005.
198. G. J. Kroes, E. J. Baerends and R. C. Mowrey, Phys. Rev. Lett., 1997, 78, 3583–3586.
199. J. Q. Dai and J. C. Light, J. Chem. Phys., 1997, 107, 1676–1679.
200. T. H. Liu, B. N. Fu and D. H. Zhang, J. Chem. Phys., 2013, 139, 184705.
201. G. J. Kroes and C. Díaz, Chem. Soc. Rev., 2016, 45, 3658–3700.
202. H. Shi, T. H. Liu, Y. Fu, X. Lu, B. Fu and D. H. Zhang, J. Phys. Chem. C, 2021, 125, 23105–23114.
203. Z. J. Zhang, T. H. Liu, B. N. Fu, X. M. Yang and D. H. Zhang, Nat. Commun., 2016, 7, 11953.
204. T. H. Liu, J. Chen, Z. J. Zhang, X. J. Shen, B. N. Fu and D. H. Zhang, J. Chem. Phys., 2018, 148, 144705.
205. X. J. Shen, Z. J. Zhang and D. H. Zhang, J. Chem. Phys., 2016, 144, 101101.
206. X. J. Shen, Z. J. Zhang and D. H. Zhang, J. Chem. Phys., 2017, 147, 024702.
207. A. K. Tiwari, S. Nave and B. Jackson, Phys. Rev. Lett., 2009, 103, 253201.
208. A. K. Tiwari, S. Nave and B. Jackson, J. Chem. Phys., 2010, 132, 134702.
209. X. J. Shen, Z. J. Zhang and D. H. Zhang, Phys. Chem. Chem. Phys., 2015, 17, 25499–25504.
210. B. Jackson and S. Nave, J. Chem. Phys., 2011, 135, 114701.
211. S. Nave and B. Jackson, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 233408.
212. H. Chadwick, H. Guo, A. Gutiérrez-González, J. P. Menzel, B. Jackson and R. D. Beck, J. Chem. Phys., 2018, 148, 014701.
213. H. Guo, J. P. Menzel and B. Jackson, J. Chem. Phys., 2018, 149, 244704.
214. H. Guo and B. Jackson, J. Chem. Phys., 2016, 144, 184709.
215. B. Jackson, F. Nattino and G. J. Kroes, J. Chem. Phys., 2014, 141, 054102.
216. V. L. Campbell, N. Chen, H. Guo, B. Jackson and A. L. Utz, J. Phys. Chem. A, 2015, 119, 12434–12441.
217. H. Guo, A. Farjamnia and B. Jackson, J. Phys. Chem. Lett., 2016, 7, 4576–4584.
218. R. Welsch, K. Song, Q. Shi, S. C. Althorpe and T. F. Miller III, J. Chem. Phys., 2016, 145, 204118.
219. S. Habershon, D. E. Manolopoulos, T. E. Markland and T. F. Miller III, Annu. Rev. Phys. Chem., 2013, 64, 387–413.
220. Q. H. Liu, L. Zhang, Y. Li and B. Jiang, J. Phys. Chem. Lett., 2019, 10, 7475–7481.
221. C. Li, Q. Liu, L. Zhang, Y. Li and B. Jiang, Mol. Phys., 2021, 120, e1941367.
222. L. Zhang, J. Zuo, Y. V. Suleimanov and H. Guo, J. Phys. Chem. Lett., 2023, 14, 7118–7125.
223. H. Y. Jiang, M. Kammler, F. Z. Ding, Y. Dorenkamp, F. R. Manby, A. M. Wodtke, T. F. Miller III, A. Kandratsenka and O. Bünermann, Science, 2019, 364, 379–382.
224. K. Gu, C. Li, S. Lin and H. Guo, J. Phys. Chem. C, 2022, 126, 17093–17101.
225. C. Li, Y. Li and B. Jiang, Chem. Sci., 2023, 14, 5087–5098.
226. M. Kurahashi and Y. Yamauchi, Phys. Rev. Lett., 2013, 110, 246102.
227. L. Vattuone, M. Rocca, C. Boragno and U. Valbusa, J. Chem. Phys., 1994, 101, 713–725.
228. A. Raukema, D. A. Butler and A. W. Kleyn, J. Phys.: Condens. Matter, 1996, 8, 2247–2263.
229. M. Kurahashi, J. Chem. Phys., 2019, 151, 084702.
230. M. Minniti, D. Farías, P. Perna and R. Miranda, J. Chem. Phys., 2012, 137, 074706.
231. D. Y. Zhang, C. Jansen, A. W. Kleyn and L. B. F. Juurlink, Phys. Chem. Chem. Phys., 2023, 25, 14862–14868.
232. J. Hall, O. Saksager and I. Chorkendorff, Chem. Phys. Lett., 1993, 216, 413–417.
233. M. Yata and H. Rouch, Appl. Phys. Lett., 1999, 75, 1021–1023.
234. P. Junell, M. Ahonen, M. Hirsimäki and M. Valden, Surf. Rev. Lett., 2004, 11, 457–461.
235. Q. Shen, J. Wu, F. Zhou, Y. Song, W. Dong, X. Wang, T. Wang and X. M. Yang, Rev. Sci. Instrum., 2022, 93, 013201.
236. S. K. Singh and P. R. Shirhatti, J. Chem. Phys., 2024, 160, 024702.
237. H. Hou, Y. Huang, S. J. Gulding, C. T. Rettner, D. J. Auerbach and A. M. Wodtke, Science, 1999, 284, 1647–1650.
238. H. Hou, C. T. Rettner, D. J. Auerbach, Y. Huang, S. J. Gulding and A. M. Wodtke, Faraday Discuss., 1999, 113, 181–199.
239. R. Cooper, I. Rahinov, C. Yuan, X. M. Yang, D. J. Auerbach and A. M. Wodtke, J. Vac. Sci. Technol., A, 2009, 27, 907–912.
240. J. Geweke, P. R. Shirhatti, I. Rahinov, C. Bartels and A. M. Wodtke, J. Chem. Phys., 2016, 145, 054709.
241. H. F. Berger and K. D. Rendulic, Surf. Sci., 1991, 253, 325–333.
242. H. F. Berger, PhD thesis, Technische Universität Graz, 1992.
243. V. Shukla, Y. Jiao, C. Frostenson and P. Hyldgaard, J. Phys.: Condens.Matter, 2022, 34, 025902.
244. V. Shukla, Y. Jiao, J.-H. Lee, E. Schröder, J. B. Neaton and P. Hyldgaard, Phys. Rev. X, 2022, 12, 041003.
245. J. W. Song, K. Yamashita and K. Hirao, J. Chem. Phys., 2011, 135, 071103.
246. D. Chakraborty, K. Berland and T. Thonhauser, J. Chem. Theory Comput., 2020, 16, 5893–5911.
247. A. Raghav, R. Maezono, K. Hongo, S. Sorella and K. Nakano, J. Chem. Theory Comput., 2023, 19, 2222–2229.
248. K. Nakano, S. Sorella, D. Alfè and A. Zen, J. Chem. Theory Comput., 2024, 20, 4591–4604.
249. S. Pathak and L. K. Wagner, J. Chem. Phys., 2018, 149, 234104.
250. J. T. Krogel and F. A. Reboredo, J. Chem. Phys., 2020, 153, 104111.
251. T. Ichibha, Y. Nikaido, M. C. Bennett, J. T. Krogel, K. Hongo, R. Maezono and F. A. Reboredo, J. Chem. Phys., 2023, 159, 164114.
252. A. Zen, J. G. Brandenburg, A. Michaelides and D. Alfè, J. Chem. Phys., 2019, 151, 134105.
253. T. A. Anderson and C. J. Umrigar, J. Chem. Phys., 2021, 154, 214110.
254. J. P. Perdew, V. N. Staroverov, J. Tao and G. E. Scuseria, Phys. Rev. A: At., Mol., Opt. Phys., 2008, 78, 052513.
255. P. Borlido, M. A. L. Marques and S. Botti, J. Chem. Theory Comput., 2018, 14, 939–947.
256. J. W. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao, A. Ruzsinszky, G. I. Csonka, G. E. Scuseria and J. P. Perdew, Phys. Rev. Lett., 2013, 111, 106401.
257. M. M. Maier, A. V. Arbuznikov and M. Kaupp, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, 9, e1378.
258. J. W. Sun, B. Xiao and A. Ruzsinszky, J. Chem. Phys., 2012, 137, 051101.
259. A. V. Krukau, G. E. Scuseria, J. P. Perdew and A. Savin, J. Chem. Phys., 2008, 129, 124103.
260. J. Kirkpatrick, B. McMorrow, D. H. P. Turban, A. L. Gaunt, J. S. Spencer, A. G. D. G. Matthews, A. Obika, L. Thiry, M. Fortunato, D. Pfau, L. Román Castellanos, S. Petersen, A. W. R. Nelson, P. Kohli, P. Mori-Sánchez, D. Hassabis and A. J. Cohen, Science, 2021, 374, 1385–1389.
261. E. G. d'Agliano, P. Kumar, W. Schaich and H. Suhl, Phys. Rev. B: Solid State, 1975, 11, 2122–2143.
262. P. Nieto, E. Pijper, D. Barredo, G. Laurent, R. A. Olsen, E. J. Baerends, G. J. Kroes and D. Farías, Science, 2006, 312, 86–89.
263. J. Behler, J. Chem. Phys., 2011, 134, 074106.
264. B. Smits and M. F. Somers, J. Chem. Phys., 2021, 154, 074710.
265. J. Meyer and K. Reuter, Angew. Chem., Int. Ed., 2014, 53, 4721–4724.
266. H. W. Sheng, M. J. Kramer, A. Cadien, T. Fujita and M. W. Chen, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 134118.
267. T. J. Frankcombe, M. A. Collins and D. H. Zhang, J. Chem. Phys., 2012, 137, 144701.
268. G. J. Kroes, J. G. Snijders and R. C. Mowrey, J. Chem. Phys., 1995, 103, 5121–5136.
269. J. Daru, H. Forbert, J. Behler and D. Marx, Phys. Rev. Lett., 2022, 129, 226001.
270. Y. Litman, E. S. Pós, C. L. Box, R. Martinazzo, R. J. Maurer and M. Rossi, J. Chem. Phys., 2022, 156, 194106.
271. Y. Litman, E. S. Pós, C. L. Box, R. Martinazzo, R. J. Maurer and M. Rossi, J. Chem. Phys., 2022, 156, 194107.
272. R. Cooper, Z. S. Li, K. Golibrzuch, C. Bartels, I. Rahinov, D. J. Auerbach and A. M. Wodtke, J. Chem. Phys., 2012, 137, 064705.
273. R. Tibshirani, J. Roy. Stat. Soc. B, 1996, 58, 267–288.
274. M. A. Efroymson, Multiple regression analysis, Mathematical methods for digital computers, Wiley, New York, 1960.
275. P. Morgante and R. Peverati, Phys. Chem. Chem. Phys., 2019, 21, 19092–19103.
276. M. Korth and S. Grimme, J. Chem. Theory Comput., 2009, 5, 993.
277. H. S. Yu, W. J. Zhang, P. Verma, X. He and D. G. Truhlar, Phys. Chem. Chem. Phys., 2015, 17, 12146–12160.
278. H. S. Yu, X. He and D. G. Truhlar, J. Chem. Theory Comput., 2016, 12, 1280–1293.
279. H. S. Yu, X. He, S. L. Li and D. G. Truhlar, Chem. Sci., 2016, 7, 5032–5051.
280. S. J. Klippenstein, V. S. Pande and D. G. Truhlar, J. Am. Chem. Soc., 2014, 136, 528–546.
281. S. Dahl, A. Logadottir, R. C. Egeberg, J. H. Larsen, I. Chorkendorff, E. Törnqvist and J. K. Nørskov, Phys. Rev. Lett., 1999, 83, 1814–1817.
282. T. Zambelli, J. Wintterlin, J. Trost and G. Ertl, Science, 1996, 273, 1688–1690.
283. K. Golibrzuch, P. R. Shirhatti, J. Geweke, J. Werdecker, A. Kandratsenka, D. J. Auerbach, A. M. Wodtke and C. Bartels, J. Am. Chem. Soc., 2015, 137, 1465–1475.
284. D. J. Harding, J. Neugebohren, H. Hahn, D. J. Auerbach, T. N. Kitsopoulos and A. M. Wodtke, J. Chem. Phys., 2017, 147, 013139.
285. J. Neugebohren, D. Borodin, H. W. Hahn, J. Altschäffel, A. Kandratsenka, D. J. Auerbach, C. T. Campbell, D. Schwarzer, D. J. Harding, A. M. Wodtke and T. N. Kitsopoulos, Nature, 2018, 558, 280–283.
286. G. B. Park, T. N. Kitsopoulos, D. Borodin, K. Golibrzuch, J. Neugebohren, D. J. Auerbach, C. T. Campbell and A. M. Wodtke, Nat. Rev. Chem., 2019, 3, 723–732.
287. D. Borodin, K. Golibrzuch, M. Schwarzer, J. Fingerhut, G. Skoulatakis, D. Schwarzer, T. Seelemann, T. N. Kitsopoulos and A. M. Wodtke, ACS Catal., 2020, 10, 14056–14066.
288. D. Borodin, I. Rahinov, O. Galparsoro, J. Fingerhut, M. Schwarzer, K. Golibrzuch, G. Skoulatakis, D. J. Auerbach, A. Kandratsenka, D. Schwarzer, T. N. Kitsopoulos and A. M. Wodtke, J. Am. Chem. Soc., 2021, 143, 18305–18316.
289. D. Borodin, N. Hertl, G. B. Park, M. Schwarzer, J. Fingerhut, Y. Q. Wang, J. X. Zuo, F. Nitz, G. Skoulatakis, A. Kandratsenka, D. J. Auerbach, D. Schwarzer, H. Guo, T. N. Kitsopoulos and A. M. Wodtke, Science, 2022, 377, 394–398.
290. M. Schwarzer, N. Hertl, F. Nitz, D. Borodin, J. Fingerhut, T. N. Kitsopoulos and A. M. Wodtke, J. Phys. Chem. C, 2022, 126, 14500–14508.

---