First-principles full-dimensional modelling of vibrational energy transfer of molecule scattering from metal surfaces

Domain: physics
Published: 2025-11-26
Editorial outcome: minor_revision

Abstract
Energy transfer during molecular collisions on metal surfaces plays a pivotal role in a host of critical interfacial processes. Despite significant efforts, our understanding of relevant energy transfer mechanisms, even in an extensively-studied benchmark like NO scattering from Au(111), remains far from complete. To fully disentangle different energy transfer channels, we develop a first-principles nonadiabatic dynamical model that incorporates explicitly all degrees of freedom and the interfacial electron transfer. Our simulations successfully reproduce most experimental observations on vibrational relaxation and excitation of NO molecules under varying initial conditions. The observed steric effect varying with the initial vibrational state is understood by the change of orientational dependence of the metal-to-molecule electron transfer. This model also identifies that translational motion could couple to molecular vibration directly, while the translation-to-electron nonadiabatic coupling is not significant. These valuable insights highlight the importance of treating both adiabatic and nonadiabatic energy transfer pathways on equal footing, offering significant implications for modeling energy transfer processes in more complex systems, such as plasmonic photocatalysis.

Introduction
Chemical reactions at metal surfaces have profound implications for interfacial applications like heterogeneous catalysis, crystal growth, and etching. Vibrational energy transfer (VET) in molecule-metal surface encounters is highly relevant to surface reactions, as molecular vibration corresponds exactly to the motion of bond breaking/forming1–4. Vibrational energy can be transferred to molecular translation, rotation, and surface phonons, an adiabatic process that can be described within the Born-Oppenheimer approximation (BOA)5,6. Beyond that, VET on metal surfaces is more complicated due to the nonadiabatic energy exchange between molecular vibration and electron-hole pairs (EHPs) of the metal, representing a breakdown of BOA7–9. There has been ample evidence of this nonadiabatic VET channel, causing vibrational excitation/deexcitation, electron emission, or short vibrational lifetimes10–15. Such nonadiabatic processes are also closely related to plasmonic catalysis, where hot electrons generated by plasmon decay may strongly couple to molecular vibration, resulting in vibrationally excited molecules that are more prone to dissociation than thermal ones16–19. An in-depth understanding of the VET dynamics at metal surfaces is thus highly desirable. However, our understanding of different VET channels remains far from complete, due largely to the lack of an accurate theoretical description of both adiabatic and nonadiabatic processes on equal footing. NO scattering from Au(111) represents one of the most important examples in this aspect, attracting extensive experimental and theoretical attention8,12,13,15,20–37. State-to-state measurements by Wodtke and coworkers have revealed a number of fascinating features, including the multi-quantum vibrational relaxation of NO(vi > 0) and its dependence on the incidence translational energy (Ei) and orientation13,25,27,29, the vibrational excitation of NO(vi = 0) and its dependence on surface temperature (Ts) and Ei15,30. Such abundant experimental data have provided detailed information on how molecular vibration couples to other degrees of freedom (DOFs), serving as valuable benchmarks to stringently test nonadiabatic theories20–22,24,31,33,34,36–40.

Despite substantial efforts in the past two decades, no theory has so far been capable of quantitatively reproducing most experimental observations for NO collisions on Au(111)1. The greatest challenge is to describe the nonadiabatic interactions between the molecular states and the metallic continuum in full dimensionality. To make the simulation tractable, most theoretical studies have relied on mixed quantum-classical approaches, which treat the nuclear motion classically while electronic motion quantum mechanically. Their coupling is either addressed perturbatively on top of the ground state, as in molecular dynamics with electronic friction (MDEF)20,33,41, or explicitly by considering two local electronic configurations as a result of electron transfer, as in the independent electron surface hopping (IESH)21–23,37 and broadened classical master equation methods36,42. Fully quantum dissipative20,40 and semi-classical39 approaches have also been applied to this system, though limited to only two nuclear DOFs. Indeed, only MDEF and IESH models have ever been integrated with realistic full-dimensional (FD) potential energy surfaces (PESs) to enable meaningful comparisons with experimental data15,21,33,37. Such comparisons based on empirical PESs parameterized by Roy, Shenvi, and Tully (referred to as RST PESs hereafter) suggested an exclusive electron transfer mediated VET mechanism, which can be characterized by IESH yet not MDEF15,23,37. However, both methods failed to capture the translational and steric dependence of vibrational inelasticity of NO(vi = 3), even qualitatively25,37,43. These failures were attributed to the inaccuracies of the RST PESs, presumably arising from the approximate charge-transfer states and less expressive empirical functions25,37.

More recently, we developed a more reliable ground state FD-PES by embedded atom neural network (EANN) fitting of thousands of density functional theory (DFT) data, on which Born-Oppenheimer molecular dynamics (BOMD) predicted very efficient adiabatic VET32,35. MDEF simulations on this PES better described the single-quantum vibrational relaxation, but still largely underestimated the extent of multi-quantum vibrational relaxation33. To go beyond the friction model, we applied a constrained DFT (CDFT) scheme44,45 to determine charge-transfer (diabatic) states of molecules at metal surfaces, which was found to more reasonably describe the electron transfer between CO and NO molecules and metal surfaces46,47. The combination of CDFT, EANN, and IESH methods enables a first-principles description for nonadiabatic effects of molecular scattering from metal surfaces. In this work, we demonstrate the success of this strategy in NO scattering from Au(111) using newly developed NN-based first-principles charge-transfer state PESs. With an appropriate density functional, IESH simulations on NNPESs achieved very good agreement with a wider range of experimental observations of NO on Au(111). Our results emphasize the necessity of properly characterizing both adiabatic and nonadiabatic energy transfer channels, towards a quantitative theoretical description for the VET dynamics of molecules at metal surfaces. The remaining discrepancies with experimental data are also discussed.

Results
Comparison of full-dimensional and two-dimensional dynamics
In Fig. 1, vibrational state distributions for scattered NO obtained by various levels of theory, for initial states of vi = 3, 11, and 16, are compared with experiments27,37. For fair comparison, unless stated elsewhere, our results are analyzed based on single-bounce trajectories, consistent with prior theoretical work and the direct scattering nature of this process as observed in experiments13. For clarity, our results based on all trajectories are shown in Supplementary Fig. 6, and no significant differences are found.

We first evaluate the accuracy of IESH by benchmarking IESH with a fully quantum method, namely hierarchical equations of motion (HEOM)40,48. Since FD HEOM simulations are currently unfeasible, we compare IESH and HEOM results based on a 2D model. As shown in Fig. 1a–c, IESH agrees well with HEOM, although the latter predicts slightly stronger vibrational inelasticity than the former in the case of NO(vi = 3). The good agreement between them largely validates the reliability of IESH.

To demonstrate the necessity of encompassing all DOFs in studying this energy transfer process, we also compare IESH results using 2D and FD models. As shown in Fig. 1a–c, the FD-IESH results are quite different from their 2D counterparts, manifesting much hotter vibrational state distributions for NO(vi = 11 and 16). Indeed, these 2D results are all based on 2D empirical PESs fitted to a limited number of CDFT points reported in ref. 46, where the NO molecule lies perpendicularly on the hcp site of Au(111) with the N atom pointing down. Electron transfer between the molecule and the surface more easily undergoes in this geometry compared to other sites and orientations (see Supplementary Fig. 7). Apparently, the effects of molecular orientation, surface site, and surface phonons are all neglected in 2D models. The large discrepancy between FD and 2D results underscores the necessity of employing FD-PESs for accurate theoretical modeling. While 2D models may coincidentally align with experimental data, our FD-IESH results remarkably underestimate the degree of vibrational relaxation, suggesting that the revPBE-based PES yields overly repulsive NO/Au(111) interactions. These findings caution against direct comparisons between experimental data and low-dimensional theoretical models, as the latter may fail to capture the true molecule-surface interactions and dynamics.

Influence of PES on VET dynamics
To study the influence of the PES on the adiabatic and nonadiabatic VET dynamics, we have constructed another set of FD-PESs based on the PW91 functional, where the ground state PES was reported previously32,33,35 and diabatic state PESs are obtained in this work as described in the “Methods” section. Figure 1d–f displays results obtained on the PW91-based PESs, where MDEF results with orbital-dependent friction (MDEF(ODF)) are taken from ref. 33. As discussed previously32,33,35, PW91-based BOMD results already predict remarkable vibrational relaxation. The MDEF(ODF) model further increases the probability of single-quantum vibrational relaxation (i.e., Δv = 1) compared to BOMD on the PW91-based PESs (Fig. 1d–f), leading to improved agreement with experimental data. However, it falls short of improving the multi-quantum vibrational relaxation (i.e., Δv > 1)33. Its failure was attributed to the invoked Markovian approximation, which neglects EHP excitations far away from the Fermi level33.

Encouragingly, the IESH model, which treats both electron transfer and EHP excitations on equal footing, achieves the best agreement with experimental observations across all initial conditions to date. Quantitatively, not only do the probability of the two-quanta relaxation (vi = 3 → vf = 1) increase remarkably, but the overall vibrational state distributions of NO(vi = 11 and 16) also shift to lower states, both of which align more closely with the experimental distributions. By comparing final vibrational state distributions predicted by BOMD and IESH, we conclude that adiabatic and nonadiabatic contributions to the vibrational relaxation with low and high initial vibrational states are comparable. In contrast, in Tully’s work21, they stated that the nonadiabatic channel dominates the VET for the NO/Au(111) scattering, which differs a lot from our conclusion, mainly resulting from their less accurate ground and charge-transfer state PESs.

In comparison, PW91-based results represent significant improvements over the revPBE-based and previous RST-based ones23,37 (Fig. 1g–i). On one hand, the adiabatic VET is much more pronounced on the PW91 PES than on the other two PESs (Note that the RST-based BOMD results exhibited minimal vibrational relaxation21, thus not shown here). This adiabatic VET arises from mode softening as the vibrating molecule approaches the dissociation barrier32,49. As depicted in Fig. 2, the barrier for NO dissociation on the adiabatic PW91 PES is 0.69 eV lower than that on the revPBE PES, while the RST PES is purely repulsive. In addition, the barrier locates “earlier” on the PW91 PES, featuring a shorter N–O distance, a higher molecule-surface distance, and a wider entrance to the barrier region (see Fig. 2b, c and Supplementary Table 1). These characteristics allow trajectories to access the barrier more readily on the PW91 PES, as demonstrated by an exemplary trajectory of NO(vi = 16) projected onto the 2D PES in Fig. 2. This facilitates significantly more efficient adiabatic VET compared to the other two PESs. On the other hand, in IESH models, the likelihood of electronic excitation is qualitatively linked to the proximity to the crossing seam between neutral and anionic states, where nonadiabatic coupling is strong and electron transfer can explicitly occur, as marked in Fig. 2 by dashed lines. Notably, trajectories on PW91 PESs are more likely to enter this strong-coupling region, making electron transfer possible and channeling vibrational energy dissipation into surface EHPs. Indeed, as compared in Supplementary Tables 2 and 3, the mean vibrational energy losses predicted by PW91-based PESs, due to both adiabatic and nonadiabatic channels, are larger than those predicted by revPBE-based ones. These results reveal that the PES landscape governs simultaneously both adiabatic and nonadiabatic VET channels. This is achieved by modulating the molecular accessibility to the dissociation barrier and state-crossings, which are typically located in similar regions.

Since the PW91-based IESH model appears to describe well both adiabatic and nonadiabatic VET processes, we next extend it to more comprehensive scenarios, including the dependence of vibrational relaxation probability on the incidence energy and initial orientation distribution, the vibration-to-translation coupling, and the vibrational excitation of NO(vi = 0). Some of these phenomena have never been adequately described by any previous theoretical models.

Incidence energy dependence of NO scattering
Figure 3 shows the experimentally measured product state branching ratios of NO(vi = 2 and 3) as a function of Ei. The observed vibrational inelasticity increases monotonically as Ei rises12,25. PW91-based BOMD results qualitatively capture this Ei-dependence, but significantly underestimate absolute ratios. MDEF(ODF) results show some improvement in the vibrationally elastic (vi = 3) and single-quantum relaxed (vf = 2) channels, but have a minimal impact on the two-quanta relaxed (vf = 1) channel. Impressively, IESH results now well match nearly all branching ratios in the entire range of Ei, demonstrating excellent agreement with experimental data. It is easily understood that the higher incidence energy allows the NO molecule to get closer to the barrier and state-crossing point, thereby amplifying both adiabatic and nonadiabatic vibrational energy losses. In contrast, previous RST-based IESH results predict a more pronounced vibrational relaxation at lower Ei and thus much weaker Ei-dependence of the vibrational inelasticity, because of the overly softness of the RST PES.

Orientational dependence of NO scattering
An intriguing yet not fully understood experimental observation is the orientational dependence of vibrational inelasticity for the scattering of NO(vi = 3), which gradually diminishes as the molecule becomes more vibrationally excited27,29. In the experiment, the vibrationally excited NO molecule can be initially manipulated to prefer either an N-first or an O-first orientational distribution before colliding on the surface, as illustrated in Fig. 4a. Figure 4b shows that the former leads to obviously stronger vibrational relaxation than the latter during NO(vi = 3) scattering. However, as the initial vibrational state increases, the orientational dependence largely weakens at NO(vi = 11) (Fig. 4c) and nearly vanishes at NO(vi = 16) (Fig. 4d). Our PW91-based IESH results properly capture this steric effect and its variation with the initial vibrational state. In Supplementary Fig. 9, we further quantify the adiabatic and nonadiabatic contributions to the steric effect by comparing BOMD and IESH results with experimental data for NO(vi = 3). We assess the degree of vibrational relaxation using the average vibrational quantum number of scattered NO molecules, denoted as <vf>, and quantify the steric effect as the difference of <vf> between the two initial orientations, namely <vf(O-first)> - <vf(N-first)>. A larger difference indicates a stronger steric effect. The BOMD prediction of this difference is ~0.13, suggesting a minor orientational dependence due to the adiabatic VET channel. While the IESH prediction increases this value to ~0.21, in close accord with the experimental value of ~0.23. These results demonstrate that both adiabatic and nonadiabatic channels contribute to the steric dependence of VET.

To elucidate the underlying mechanism, we first check the polar angle distributions at the impact point at which the initially N-first and O-first oriented molecules turn around (see Supplementary Fig. 10). While some O-first oriented NO molecules reorient towards the N-first orientation, the polar angle distributions in the two cases remain distinct, which is the prerequisite for the stereodynamics. However, these distributions are insensitive to the initial vibrational state, failing to explain the reduced steric effect at higher vibrational states. We then resort to searching the crossings between two diabatic state PESs, plotted as a function of ZNO in Supplementary Fig. 10, with the polar angle θ fixed at 68° and 112°—values that correspond to the mean angles (<θ>) for the N-first and O-first distributions, respectively. At rNO = 1.32 Å, approximately the outer classical turning point of NO(vi = 3), the crossing point of two diabatic states shifts to higher energies as θ increases. This indicates that N-first oriented molecules would more likely undergo charge transfer than O-first ones and experience more efficient nonadiabatic VET, consistent with their stronger vibrational inelasticity upon collisions. By contrast, when rNO elongates to the turning point of NO(vi = 16), i.e., 1.59 Å, the crossings of two diabatic states become very similar in both orientations. This is in accord with the observed much weaker steric effects for high initial vibrational states. It is worth noting that previous RST-based IESH results turn out to predict no steric effect at all for NO(vi = 3) scattering, due to too strong dynamical steering21,43,50 (see Fig. 4 of ref. 21 and Supplementary Fig. 10 in the Supplementary Information, which presents similar polar angle distributions at the impact point from different initial orientations).

Vibration-to-translation coupling
Dynamical results discussed have so far primarily focused on the energy exchange between molecular vibration and surface EHPs (V–EHP). However, the energy transfer process is more complex, involving other DOFs, like molecular translation and rotation. Experimental evidence suggests that the translational motion is not merely a spectator during the electronically nonadiabatic VET process, but its precise role remains elusive. For example, as shown in Fig. 5a, the mean final translational energy (<ET>) is incrementally higher in the vibrationally de-excited channels than the elastic channel26,28. Experimentalists have proposed two possible mechanisms: a direct mechanical coupling between vibration and translation (V–T), and/or an EHP-mediated coupling, where molecular vibration and translation do not directly couple to each other but both couple to EHPs (V–EHP and EHP–T)26,28. However, no theory has yet been able to determine which mechanism dominates.

Figure 5 compares the calculated and experimental <ET> of scattered NO molecules, as a function of final rotational energy (ER) for three product channels (vi = 3 → vf = 3, 2, 1), at Ei = 0.98 eV. The experimentally observed <ET> decreases linearly with increasing ER for all three channels, exhibiting an anticorrelation between outbound translation and rotational excitation (T–R). In addition, the evidence of V–T energy transfer is that <ET,vf=3> extrapolated to ER = 0 is lower by ~0.08 eV than <ET,vf=2> and ~0.12 eV than <ET,vf=1>, smaller than the corresponding vibrational excitation energy, respectively. In comparison, PW91-based IESH results not only correctly predict the anticorrelation between <ET> and ER, but also reproduce reasonably well the amount of V–T energy transfer, e.g., ~0.09 eV for <ET,vf=2>-<ET,vf=3> and ~0.15 eV for <ET,vf=1>-<ET,vf=3>, although the slope of d<ET>/dER is more negative, implying the anisotropy of the entrance of the PES is less well described.

Interestingly, Fig. 5c shows that PW91-based BOMD results also exhibit a qualitatively similar anticorrelation between <ET> and ER across three well-separated vibrational channels, suggesting a direct V–T coupling with no need to invoke nonadiabatic effects. Indeed, the differences in <ET> between adjacent vibrational channels are roughly constant at ~0.22 eV, corresponding approximately to the energy of a single vibrational quantum (Δv = 1). This means that vibrational energy can exclusively convert to the outbound translation because of vibrational deexcitation within BOA. When the nonadiabatic channel is open, part of the vibrational energy flows to EHPs, thereby reducing the differences of <ET> between vibrationally relaxed and elastic channels. This analysis also conforms to the mean energy changes in different DOFs during NO(vi = 3) scattering (Supplementary Table 3). Compared to BOMD, the average translational energy loss is barely changed when electronic DOFs are included via IESH, showing no evidence of energy exchange between translational and electronic DOFs. This further confirms the negligible coupling between translation and EHPs. This feature can be understood by analyzing the direction of the nonadiabatic coupling vector whenever surface hopping occurs in IESH simulations. As shown in Fig. 5d, the angular distribution between the nonadiabatic coupling vector and N–O bond (α) is sharply peaked at parallel directions, either 0° or 180°. In such scenarios, the nonadiabatic coupling vector is either parallel or antiparallel to the N–O vibrational motion, resulting in the predominance of V–EHP coupling. These results clearly demonstrate that molecular vibrational energy can transfer to both molecular translation and metallic EHPs during scattering, yet the direct T–EHP coupling is negligible, at least in the conditions considered here.

Vibrational excitation of NO(vi = 0) scattering
Apart from vibrational relaxation of NO at various initial conditions, vibrational excitation of NO(vi = 0) was also observed during its scattering from Au(111)15,30. In this process, the electronic energy is in turn transferred to molecular vibration. Indeed, this process is relevant to plasmonic catalysis, where nonadiabatic energy transfer from hot electrons into vibrationally excited states facilitates the bond cleavage19. As shown in Fig. 6, our PW91-based IESH results compare well with experimental data on the vibrational excitation probabilities for NO(vi = 0 → vf = 1), in a range of Ts from 473 to 973 K at two representative Ei values. Both results are below the thermal limit, indicating that vibrational excitation occurs in a direct scattering process. Experimental vibrational excitation probabilities can be fitted to an Arrhenius equation on Ts. PW91-based IESH results generally follow this Arrhenius relation, except for the deviation at Ei = 0.63 eV, Ts = 473 K. The overall agreement observed here is comparable to that of the earlier RST-based IESH model15.

Methods
CDFT calculations
A key difference between the present and a previous model23 is that U0 and U1 are now determined by CDFT energies with a constrained Bader charge of 0 e− and −1 e− onto the NO moiety, respectively. Their couplings are then derived by their energies and the ground state DFT energy (Eg) via,
Vc = sqrt((U0 − Eg) × (U1 − Eg)), (1)
which is then used to evaluate Vak. According to Wu and Van Voorhis44, in CDFT, the energy functional is defined by modifying the traditional DFT functional, E[ρ], with a Lagrange multiplier (V) dependent on the charge constraint,
W[ρ, V] = E[ρ] + V (integral wc(r)ρ(r)dr − Nc). (2)
Here, wc(r) is a spatial weight function of the charge density ρ(r), and Nc is the target value of the charge constraint. The energy of the charge-constrained state is obtained by self-consistently minimizing W with respect to ρ and maximizing W with respect to V,
ECDFT = minρ maxV [W[ρ, V]]. (3)
The CDFT diabatic state energies and forces were obtained using the CDFT module implemented in CP2K45—a general way to calculate charge-constrained states for molecule-metal systems46. Compared to the electric field method23, diabatic states defined by CDFT are well-behaved and smooth in the entire dynamically relevant configuration space, enabling constructions of FD diabatic PES with accuracy at the first-principles level. For the current NO/Au(111) system, the convergence criterion for the charge constraint was set to 5 × 10−3 e−. In detail, Kohn−Sham orbitals were expanded into a mixed Gaussian plane-wave basis, and the ionic cores were approximated by the Goedecker–Teter–Hutter (GTH) pseudopotentials62. Valence electrons were expanded with the DZVP-MOLOPT-SR-GTH basis set for all elements. Since NO is an open-shell molecule, all calculations were spin polarized. The surface was modeled by a 6×6 supercell with four metal layers. The top two layers of Au atoms were movable, and a vacuum space of 25 Å in the Z direction was used to separate the periodic slabs. A Fermi-Dirac smearing method was used with an electronic temperature of 1200 K to maintain the fractional occupation of states near the Fermi energy and to speed up self-consistent field convergence. Most calculations in this work have been done using the revPBE functional63, except in some test calculations discussed below, where the PW91 functional64 was used.

Ground and charge-transfer state PESs
EANN method. A well-established EANN method65,66 was employed to learn CDFT energies and forces to yield charge-transfer state PESs, in the same way as representing the ground state PES. This allows us to construct the FD-IESH Hamiltonian without introducing any empirical parameters. In this approach, the total energy of the system is decomposed as the sum of atomic energies, each of which is represented by an atomic neural network,
E = sum_i E_i = sum_i NN(ρi), (4)
where ρi are embedded atom density (EAD) features regarded as the input to each atomic neural network that can be simply expressed by the square of the linear combination of Gaussian-type orbitals (GTOs) located at neighbor atoms,
ρi = sum_{lx,ly,lz=L} sum_{j≠i} [sum_{c} c_{Z_j} Z_j φ_{α,rs}^{lx,ly,lz}(rij) fc(rij)]^2. (5)
In Eq. (5), Ncut is the total number of neighbor atoms within a cutoff radius (rc), c_{Z_j} is the element-dependent orbital coefficient which can be optimized together with NNs parameters, and fc(rij) is a cutoff function to ensure that the contribution of each neighbor atom decays smoothly to zero at rc. The GTO, φ_{α,rs}^{lx,ly,lz}(rij), is defined by,
φ_{α,rs}^{lx,ly,lz}(rij) = x^lx y^ly z^lz exp(−α|rij − rs|^2), (6)
where r = (x, y, z) is the Cartesian coordinates of the embedded atom i relative to a neighbor atom j, rij is its length, α and rs are parameters that determine radial distributions of GTO, lx + ly + lz = L specifies the orbital angular momentum (L).

The revPBE-based PESs. It should be noted that we have constructed PESs based on two density functionals to investigate the influence of the energy landscape on the adiabatic and nonadiabatic VET dynamics. We first chose the revPBE functional63 for two main reasons. First, IESH and HEOM results based on 2D (only includes the N–O bond length (rNO) and the molecular height on the surface (ZNO), as defined in Supplementary Fig. 1), PESs fitted from dozens of data points using the revPBE functional were available40. This allowed us to conveniently compare these results with our FD-IESH results to validate the accuracy of IESH and reveal the importance of multidimensional effects. Second, the revPBE functional predicts a binding energy of NO with Au(111) that agrees with the experimental value67 (see Supplementary Table 1). However, as discussed above, the revPBE-based PESs appear too repulsive to enable efficient VET. Data points used to construct the revPBE-based ground, neutral, and anionic PESs were generated by two procedures. First, approximately 500 configurations for the NO/Au(111) system were randomly generated and computed to construct an initial ground state PES. The rest 2000 configurations were iteratively selected from quasi-classical trajectories on the iteratively refined ground state PES via both geometric and energetic criteria until the dynamical results converged. After excluding the non-converged configurations in DFT calculations, a total of 2285, 2165, 2047 data points were used to fit the ground, neutral, and anionic state EANN PESs, respectively, among which 90% of data points were used for training, while the remaining 10% were reserved for validation. These EANN hyperparameters take values of L = 0–2, rc = 7.6 Å, α = 0.53 Å−2, and rs = 0.0–7.4 Å every 0.6 Å, yielding 39 EAD features. The atomic neural network comprises two hidden layers with 30 and 40 neurons for the ground, neutral, and anionic PESs. The overall RMSEs of the ground, neutral, and anionic PESs are 11.8, 59.9, and 75.8 meV for energies, and 17.7, 60.4, and 113.0 meV/Å for atomic forces, respectively. The significantly larger errors in the neutral and anionic state PESs compared to the ground state PES primarily stem from the broader energy range of diabatic states and extra errors in CDFT calculations.

The PW91-based PESs. While the potential energy landscape has been shown to critically regulate adiabatic energy transfer dynamics32,49, its role in nonadiabatic energy transfer dynamics during molecule-metal scattering remains relatively unclear. Herein, we further construct diabatic state PESs using the PW91 functional, based on which BOMD predicted a large amount of vibrational energy loss after NO scattering from Au(111)32,35. As illustrated in Supplementary Note 1, the energy differences between the neutral and ground states (U0 − Eg) and between the anionic and ground states (U1 − Eg) calculated using the revPBE and PW91 functionals are very similar. To avoid performing costly CDFT calculations for thousands of configurations required for PES fitting, we add (URE − ERE) calculated on the revPBE PES to the PW91 ground state (EPW) to construct the PW91 neutral and anionic state PESs,
UPW = EPW + (URE − ERE), (7)
0 0 g 0 g
The PW91 anionic state PES can also be obtained by,
UPW = EPW + (URE − ERE), (8)
1 1 g 1 g
where (URE − ERE) is the energy difference calculated on the revPBE PES. The PW91-based PESs were subsequently applied in IESH simulations and compared with a variety of experiments.

IESH dynamics
In this work, nonadiabatic dynamics simulations for NO scattering from Au(111) are carried out by the IESH method22,68. The method is based on a discretized version of the Newns–Anderson (NA) Hamiltonian69, which describes an impurity level (here the lowest unoccupied molecular orbital (LUMO) of the molecule) interacting with a manifold of electronic states with the assumption that the metallic electrons are non-interacting. The discretized many-electron NA Hamiltonian is written as,
Hel(R) = U0(R) + (U1(R) − U0(R))a^+a + sum_{k=1}^M ε_k c_k^+c_k + sum_{k=1}^M Vak(R)(a^+c_k + c_k^+a). (9)
Here, U0(R) is the PES for the interaction of a neutral molecule (with the LUMO orbital, |a〉, unfilled) with the metal surface, U1(R) is the PES for the interaction of a negatively charged molecule (with |a〉 filled) with the metal surface, a+(a) is the creation (annihilation) operator for the molecular orbital, c_k+(c_k) is the creation (annihilation) operator for the kth metal electron. The last term couples the molecular orbital and the metal orbital.

Due to the assumption of independent electrons, this Hamiltonian can be written as the sum of one-electron Hamiltonians, each of which is given by,
H1el(R) = (U1(R) − U0(R))|a〉⟨a| + sum_{k=1}^M ε_k |k〉⟨k| + sum_{k=1}^M Vak(R)(|a〉⟨k| + |k〉⟨a|). (10)
By diagonalizing Eq. (10), a series of eigenvalues (E) and eigenstates |φj(R)〉 are obtained. Thus, the adiabatic Hamiltonian for evolving the nuclear subsystem can be expressed with the eigenvalues,
H(R, P) = sum_{i=1}^{3N} P_i^2 / 2M + U0(R) + sum_{j=1}^{Ne} E_j(R) s_j(t), (11)
where the total energy of electrons is represented by the sum of a set of occupied single electronic states determined by a time-dependent vector s(t). The corresponding Ne-electron eigenstate of the many-electron Hamiltonian can then be obtained by populating Ne of these one-electron orbitals in a Slater determinant,
|j〉 = |φj1, φj2, ..., φjNe〉: (12)
In this representation, the ground state corresponds to filling Ne (half of the number of metal orbitals, M) electrons from the lowest single-electron orbital up to the Fermi level, while the excited states correspond to the determinants in which one or more electrons jump to empty orbitals above the Fermi level.

In IESH, the total many-electron wavefunction is a single Slater determinant state,
|Ψ〉 = |ψ1, ψ2, ..., ψNe〉: (13)
Each electron in the system is evolved independently according to the one-electron Hamiltonian H1el,
iℏ ∂|ψ_α〉 / ∂t = H1el(R)|ψ_α〉, (14)
where α goes from 1 to Ne. We expand each single-electron wave function in terms of (M+1) electronic basis functions,
ψ_α(r, R, t) = sum_j c_j^α(t) φ_j(r, R), (15)
where {c_j^α} are the complex expansion coefficients for the αth electron. Substituting Eq. (15) into Eq. (14), we get,
i c_k^α = E_k(R) c_k^α − i sum_{j≠k} R_dot d_kj(R) c_j^α, (16)
where d_kj is the nonadiabatic coupling between orbitals k and j. In practice, Eqs. (11) and (16) describe the evolution of the nuclear and electronic subsystems. Within the IESH approach, the hopping probability between electronic states represented by |k> and |j> is non-zero only when |k> and |j> differ in exactly one occupied orbital, that is, k1=j1, k2=j2, ..., kNe=jNe, except for a single orbital at which they differ (k≠j). And the probability of a hop from the occupied orbital k to an unoccupied orbital j is calculated by,
g_{k→j} = max(0, −2 Re(A_{kj} R_dot d_{jk}) Δt / A_{kk}), (17)
where Δt is the electronic time step, A_{kj} and A_{kk} are elements of the electronic density matrix. After these modifications, the IESH method follows the conventional fewest-switch surface hopping algorithm22,70.

Quasi-classical trajectory implementation
To compare with scattering experiments, we have performed quasi-classical trajectory calculations for the nuclear dynamics. The Jacobi coordinate frame of the NO/Au(111) system for describing the state-to-state scattering process is illustrated in Supplementary Fig. 1. The BOMD and IESH simulations were performed using an in-house modified version of the VENUS code71,72. In this work, all trajectories started with NO molecules positioned with the center of mass 8.0 Å above the topmost layer of Au(111), i.e., ZNO = 8.0 Å (see Supplementary Fig. 1). The lateral coordinates of the center of mass of NO (X and Y) were randomly sampled in the simulation supercell. The initial polar (θ) and azimuthal (φ) angles of the molecular internuclear vector (defined as pointing from N to O) were randomly sampled unless stated otherwise. The initial values of the internuclear distance, rNO, and the conjugated momentum were sampled semi-classically for given vibrational and rotational quantum numbers v and j73. The initial condition of surface atoms was chosen from the Anderson thermostat, but no thermostat was imposed during the collisional process. Nuclear coordinates and momenta were propagated using the Verlet algorithm with a time step of 0.1 fs, with a maximum simulation time of 10 ps. The trajectory was terminated and assigned to a “scattered” event when molecule-surface vertical distance exceeded 8.0 Å and the molecule velocity pointed away from the surface, for which the final vibrational state (vf) of NO was determined by the Einstein−Brillouin−Keller quantization74 and rotational quantum number (jf) was obtained from the quantum mechanical expression for rotational angular momentum. Alternatively, the trajectory was labeled as “trapped” if the trajectory exceeded the maximum simulation time. To obtain the final rovibrational state distributions, the fractional vibrational and rotational quantum numbers were binned into the nearest integers via a histogram binning procedure.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this paper the authors present new results obtained for the NO + Au(111) system with a focus on describing the multi-quantum vibrational energy transfer [VET] that occurs when the molecule encounters the surface. They present results using various levels of theory [wrt dynamics] including Born-Oppenheimer Molecular Dynamics [BOMD], Molecular Dynamics including Electronic Friction [MDEF] and an Independent Electron Surface Hopping [IESH] approach, using Constrained DFT [CDFT]. In the work they present results on two different multi dimensional Neural Network [NN] PESs based on a GGA level of DFT. One is based on the revPBE DFT GGA and the other on the PW91 DFT functional. The new and exciting results presented show that especially the use of the CDFT based IESH MD using the PW91 NN PES is very successful in retrieving the experimentally observed VET. The other methods, and or use of DFT functional, are in fact inferior compared to the newly found and presented key combination.\n\nAlthough a great improvement is seen in using the IESH based on CDFT on the PW91 NN PES, not all experimental observables are now correctly predicted. The authors argue that this is most likely due to issues still remaining on the GGA level [the PW91] and or perhaps the approximations used in the dynamical methods [the IESH and CDFT]. However, another possible source of error is not mentioned at all: the NN fitted PESs. At least a mention should be made of the level of accuracy and how these NN PESs have been 'trained'. NN PESs can have artifacts and issues in areas and with observables they were not well 'trained' upon and especially this is here a relevant point to address if the authors themselves also state that the success of their best results depend largely on PES details [and used GGA functional] in their discussion.\n\nMoreover, it is in my opinion that this newly found key combination [what level of dynamical theory is to be combined with what GGA DFT functional] is both a strong and a weak point at the same time. It is in fact true that the newly found 'golden combination' outperforms the others investigated, and indeed does describe the VET for the first time in an unprecedented level of accuracy, but, at the same time, it is found that it also is very depended on the ground state PES [and the DFT GGA] details. The results therefore also represent a weakness in that they are inherently very system specific, and far from generally \"ab initio\" applicable without system specific knowledge a priory. It still even remains a question if this 'golden combination' is actually the best there is, given the fact only two DFT GGA functionals have been tried. It could also still all be 'a fortuitous case' as far as is discussed in the paper now [regarding the sensitivity on the ground state PESs details].\n\nSo, having found this 'golden combination' for NO on Au(111) is by itself important [indeed the research is also also state-of-the-art] and deserves a publication in decent journal, it is however not a combination that includes a newly developed [general] method, neither is it a finding of new physics as such for this specific system or a class of systems. I recommend therefore submitting this paper to another journal where then also the room of a more in-depth discussion is possible regarding the sources of errors that do remain.", "reviewer_2": "This remarkable article addresses a long standing and central issue in the dynamics of molecular interactions at surfaces – what is the relative importance of effects that are within the Born Oppenheimer Approximation, i.e., electronically adiabatic effects and those that involve electronically nonadiabatic transitions. They address this question by using state of the art theoretical methods to do calculations modeling the rich trove of experimental work on NO vibrationally inelastic scattering from metal surfaces and compare models that that include or exclude nonadiabatic transitions but otherwise are on an equal footing with respect to the potential energy surface and dimensionality of the problem. The agreement of the best models to a broad range of experimental data as stunning and represents a tremendous advance. Particularly remarkable is the ability of this work to reproduce the subtle results on the initially molecular orientation dependance of vibrational energy transfer and its relation to crossing of the seam between different states on the potential energy surface. The authors clearly demonstrate what they state in the title: both adiabatic and nonadiabatic effects matter.\n\nThese findings represent a dramatic advance in our understanding and ability to calculate vibrationally energy transfer including nonadiabatic electronic effects that should be of broad interest.\n\nI recommend publication in Nature Communications after the authors address the remarks given in the remainder of this review.\n\n1. Abstract\nThe abstract needs work.\nIt is OK up to the 4th sentence where the authors write “Our simulations reproduce, for the first time, most experimental observations on vibrational relaxation and excitation of NO molecules under varying initial conditions and clearly elaborate the respective adiabatic and nonadiabatic contributions”.\nThis is followed by “This model identifies the direct adiabatic vibration-to-translation coupling as the dominant role of translation, while excluding significant translation-to-electron nonadiabatic coupling”.\nTaken in the context of the previous sentence the reader will likely think that that this sentence is providing elaboration on the respective adiabatic and nonadiabatic contributions. Thus the sentence could be understood as saying that direct vibration-to-translation coupling is the dominant mechanism in vibrational relaxation. I think with this sentence the authors are referring to the results on the transfer of vibrational energy to translation in vibrational deexcitation, which they show arises in an adiabatic model with no participation of ehp coupling to translational motion of the NO molecule. But this is not conveyed by saying “the model identifies the direct adiabatic vibration-to-translation coupling as the dominant role of translation.” Translation has many roles including pushing the NO to the seam in the PES, exciting phonons in the collision, etc.\nI suggest the authors add a sentence or two after “and clearly elaborate the respective adiabatic and nonadiabatic contributions” to give some elaboration. For example they might say something like vibrational relaxation of highly vibrationally excited NO has comparable contributions from adiabatic and nonadiabatic effects. Then they could say the coupling of the energy lost in vibrational deexcitation is predominantly adiabatic and is not significantly affected by nonadiabatic translation-to-electron coupling.\n\n2. Comparison of full-dimensional and two-dimensional dynamics\nThis section would benefit from some introductory text to help the reader understand the motivations and questions that are being asked. There seem to be two motivations: (1) to see how important multidimensional effects are and (2) to used 2D calculations to validate IESH since a full quantum method is possible in 2D but not currently feasible in FD. Starting the section with some text to guide the reader rather the starting with introducing a figure would help. It might be helpful to explain why was the revPBE functional chosen for this work\n\n3. Influence of PES on VET dynamics\nThis section could benefit from clearer descriptions. For example I can not understand what is meant by “while the MDEF(ODF) model further increases the probability of single-quantum vibrational relaxation (i.e. ∆v=1)”. Please refer to which panel is referred to.\nThe statements about IESH get a little lost by being tacked on at the end of a paragraph focused on other things. The agreement with IESH should be emphasized more, for example by starting a new paragraph saying at the outset it gives the best agreement with experiment to date.", "reviewer_3": "The manuscript addresses the study of a rather complex system, which remains the subject of ongoing debate: the reaction of NO with Au(111). It employs a detailed methodology aimed at accurately reproducing the involved physical phenomena, and the analysis is thorough and significant. However, the title is unclear, as it implies that adiabatic and non-adiabatic channels can be treated independently, which seems incorrect. These channels cannot be fully separated; calculations can be entirely adiabatic, but what does it mean to be completely non-adiabatic? Furthermore, while this work emphasizes the importance of non-adiabatic effects, this issue was already raised in Tully and Co.'s seminal article [Science 326, 5954 (2009)].\n\nThe comparison between the new theoretical results and experimental data is outstanding and underscores the importance of performing dynamics on a highly accurate potential energy surface. Nonetheless, there is some debate regarding the statement that the agreement cannot be due to error compensation, especially considering that the potential energy surface was calculated using the PW91 functional.\n\nOverall, the results are compelling and warrant publication. The methodology includes all necessary components to reproduce the phenomena, and the findings are meaningful. However, the conclusions may overstate the significance of the work. It might be more appropriate for publication in a specialized journal.\n\nOn the other hand, the impact of this research could be enhanced by demonstrating that the same methodology can be applied to other complex systems, such as HCl/(Au(111)). This would showcase the broader applicability and relevance of the approach."}, "author_response": "Response to Reviewers\nNCOMMS-25-27914: First-principles understanding of vibrational energy\ntransfer in molecule-surface scattering: both adiabatic and nonadiabatic\nchannels matter\nReviewer Comments: Blue\nAuthor Responses: Black\nRevisions: Red\nComments from Reviewer #1 and our detailed responses\nGeneral comments: In this paper the authors present new results obtained for the\nNO + Au(111) system with a focus on describing the multi-quantum vibrational energy\ntransfer [VET] that occurs when the molecule encounters the surface. They present\nresults using various levels of theory [wrt dynamics] including Born-Oppenheimer\nMolecular Dynamics [BOMD], Molecular Dynamics including Electronic Friction\n[MDEF] and an Independent Electron Surface Hopping [IESH] approach, using\nConstrained DFT [CDFT]. In the work they present results on two different multi\ndimensional Neural Network [NN] PESs based on a GGA level of DFT. One is based\non the revPBE DFT GGA and the other on the PW91 DFT functional. The new and\nexciting results presented show that especially the use of the CDFT based IESH MD\nusing the PW91 NN PES is very successful in retrieving the experimentally observed\nVET. The other methods, and or use of DFT functional, are in fact inferior compared to\nthe newly found and presented key combination.\nAuthor’s reply: We thank Reviewer #1 for carefully reading our manuscript and\nfor the generally positive comments. In response, we have provided detailed\nexplanations and additional results to address the remaining concerns.\nComment 1: Although a great improvement is seen in using the IESH based on\nCDFT on the PW91 NN PES, not all experimental observables are now correctly\npredicted. The authors argue that this is most likely due to issues still remaining on the\nGGA level [the PW91] and or perhaps the approximations used in the dynamical\nmethods [the IESH and CDFT]. However, another possible source of error is not\nmentioned at all: the NN fitted PESs. At least a mention should be made of the level of\naccuracy and how these NN PESs have been 'trained'. NN PESs can have artifacts and\nissues in areas and with observables they were not well 'trained' upon and especially\n1\nthis is here a relevant point to address if the authors themselves also state that the\nsuccess of their best results depend largely on PES details [and used GGA functional]\nin their discussion.”\nAuthor’s reply: We thank the reviewer for pointing out the possible error arising\nfrom the NN PESs. The level of accuracy of the PESs and the details of data selection\nand training have been described in section I.B of the Supplementary Information (SI).\nFor the revPBE-based PESs, for example, the root-mean-square-errors (RMSEs) of the\nground, neutral, and anionic PESs are 0.08, 0.41, and 0.52 meV/atom, respectively.\nThese RMSEs are much lower than typical errors of atomistic NN PESs for high-\ndimensional systems (several meV/atom), suggesting that our NN PESs have been well\ntrained. We do note that the RMSEs for the neutral and anionic state PESs are larger\nthan that of the ground state PES, which is likely due to the much broader energy range\nof the diabatic states and extra convergence difficulties in CDFT calculations.\nTo further validate the accuracy of the trained PESs, we have performed additional\ntests by randomly excluding 10% of data points (~200 points) from the original dataset\nand retraining the PESs. The test RMSEs of these retrained PESs for these removed\ndata points were 0.15, 0.56, and 0.57 meV/atom for the ground, neutral, and anionic\nstates, respectively, which are comparable to the errors of the original PESs. IESH\nsimulations based on these retrained PESs yielded very similar results to those based\non the original PESs (see Fig. R1), indicating that the error from the PES fitting is not\nsignificant.\nChanges in revised manuscript:\n(i) Fig. R1 was added in the revised Supplementary Information (SI) as Fig. S3.\nThe convergence test of PESs has been added on page S14 in the revised SI:\n“To validate the convergence of the obtained PESs, we performed additional tests\nby randomly excluding 10% of the data points from the initial dataset and retraining the\nPESs using the remaining 90%. The test RMSEs of these 10％ data points were 0.15,\n0.56, and 0.57 meV/atom for the ground, neutral, and anionic PESs, respectively—\nvalues that are comparable to the total RMSEs of the original PESs. We then carried\nout IESH simulations based on this new set of PESs. As shown in Fig. S3, the IESH\nresults using the retrained PESs are very similar to those based on the original PESs,\nshowing the satisfying convergence of our used PESs with respect to the number of\ndata points.”\n(ii) We have also added a brief error analysis related to the PES training on page\n22 of the revised manuscript:\n2\n“The PES training may also introduce additional errors to our model. However,\ngiven the small root-mean-square-errors and the good convergence behavior of the\nPESs (see sections I.B and II.B of the SI), we expect a minor influence of PES training\nerrors on the simulation results.”\nFig. R1. Convergence tests of PESs. (a) Comparison of IESH predicted final\nvibrational state distributions of NO(v = 16, E = 0.52 eV) scattering from Au(111)\ni i\nusing revPBE PESs trained from the full dataset and the reduced dataset. (b) Similar to\n(a), but using PW91 PESs.\nComment 2: Moreover, it is in my opinion that this newly found key combination\n[what level of dynamical theory is to be combined with what GGA DFT functional]\nis both a strong and a weak point at the same time. It is in fact true that the newly found\n'golden combination' outperforms the others investigated, and indeed does describe the\nVET for the first time in an unprecedented level of accuracy, but, at the same time, it is\nfound that it also is very depended on the ground state PES [and the DFT GGA] details.\nThe results therefore also represent a weakness in that they are inherently very system\nspecific, and far from generally \"ab initio\" applicable without system specific\nknowledge a priory.\nAuthor’s reply: We appreciate the reviewer’s insightful comment. However, from\na theoretical perspective, we think any detailed reaction dynamics would be dependent\non the landscape of the PES (and thus the density functional), regardless of the system.\nAlthough we have used the most advanced dynamical theory currently available with\nan affordable cost, the corresponding results unavoidably depend on the PES or the\ndensity functional. This functional dependence is not a limitation specific to our study,\nbut a general phenomenon in chemical dynamics at surfaces. Even in much simpler\nsystems such as H /Cu(111), where phononic and electronic excitations play only a\n2\nminor role, the H dissociation and nonreactive scattering probabilities calculated on\n2\n3\nthe RPBE and PW91 PESs differ significantly, underscoring the functional dependence\nof scattering dynamics [Science 326, 832 (2009)]. To solve this issue, a so-called\nspecific reaction parameter (SRP) density functional—SRP48—was found to be\nsufficiently accurate in reproducing most experimental observables, which was\nhowever obtained by empirically scaling theoretical results to match experimental\nsticking probabilities and was not transferable. Similar sensitivity has been reported in\nother gas-surface systems, including HCl/Au(111) [e.g., J. Chem. Phys. 146, 164706\n(2017); J. Phys. Chem. C 124, 15944 (2020)] and H O/Ni(111) [e.g., Science 344, 504\n2\n(2014); Phys. Chem. Chem. Phys. 18, 21817 (2016)], where different functionals\nyielded different dynamical results, and no SRP functional was attempted. In short, to\nthe best of our knowledge, for a given system, at least so far, no combination of methods\ncan be claimed as generally ‘ab initio’ without prior knowledge about the system or\nwithout performing density functional tests.\nMore importantly, aforementioned systems involve no significant nonadiabatic\neffects, which have thus been often investigated on a single ground state PES. It is\nrelatively easy to elaborate how the adiabatic dynamics depends on the PES or density\nfunctional. In contrast, our study goes way beyond this by studying a system with strong\ncharacteristics of electron transfer using a state-of-the-art full-dimensional nonadiabatic\nmodel that combines CDFT, EANN, and IESH, minimizing empirical parameters in the\nPESs or dynamical models typically involved in such simulations. This enables us, for\nthe first time in this community, to directly examine how different density functionals\n(and the resulting PES landscapes) affect the strongly nonadiabatic energy transfer\nprocesses. Additionally, it should be emphasized that this ‘golden combination’ referred\nby the referee was naturally obtained by the fact that the PW91 functional describes\nwell the adiabatic and nonadiabatic energy transfer channels, rather than an artificial\nresult empirically scaled to match the experimental data. Indeed, adjusting empirical\nmodels may improve the agreement with experiments in some circumstances, but is\nnearly infeasible to reproduce a diverse range of experiments, as has been achieved in\nthis work.\nOverall, the dependence of dynamical results on the density functional should not\nbe regarded as a weak point of our study. Rather, we chose not only to only present\nthese good results in agreement with experiments, but also to honestly show a minor\nset of unsatisfactory results, discuss their relations to the PESs and suggest possible\ndirections for further improving the description of this challenging system. It is our\nopinion that these discussions make this work more comprehensive and more beneficial\n4\nto the community.\nChanges in revised manuscript: We accordingly added some discussions on page\n22 in the revised manuscript:\n“…… A combination of these state-of-the-art techniques minimizes empirical\nparameters in the PESs and dynamical models typically involved in nonadiabatic\nsimulations. This enables us, for the first time, to explicitly examine how the potential\nenergy landscape simultaneously affects adiabatic and electron transfer mediated\nnonadiabatic energy transfer processes. Through systematic comparisons of previous\n2D models and current FD models based on two different functionals, we find that both\nadiabatic and nonadiabatic energy transfer dynamics are dependent on the potential\nenergy landscape and model dimensionality, which are crucial for quantitatively\nexplaining experimental observations.”\nComment 3: It still even remains a question if this 'golden combination' is actually\nthe best there is, given the fact only two DFT GGA functionals have been tried. It could\nalso still all be 'a fortuitous case' as far as is discussed in the paper now [regarding the\nsensitivity on the ground state PESs details].\nAuthor’s reply: We understand the reviewer’s concern, but we do not claim that\nour combination of methods is the best. Neither is our aim to identify the optimal\ndensity functional for the NO/Au(111) system. The fact is that our PW91-based model\ndoes yield significantly improved results under various conditions over all previous\ntheoretical models and explain several experimental phenomena that had never been\naddressed by theory, such as the origin of vibrational-to-translational coupling and the\ndependence of steric effects on the initial vibrational state. This achievement is no doubt\na breakthrough after many years’ theoretical efforts. Of course, there remains room for\nfurther improvement of the model, as discussed in the Discussion section of the\nmanuscript.\nThe use of the two GGA density functionals are not arbitrary but well justified.\nOne the one hand, revPBE represents a more repulsive interaction, characterized by a\nshallow adsorption well (-0.14 eV vs -0.24 eV measured experimentally [Surf. Sci. 650,\n11 (2016)]) and a high dissociation barrier (3.58 eV) for NO on Au(111). Additionally,\nthis functional was previously used to fit the 2D PESs that allow us to compare the\nIESH method with others and highlight the importance of model dimensionality. On the\nother hand, PW91 yields a more attractive interaction, indicated by a deeper well (-0.36\neV) and a lower barrier (2.88 eV) (see Table S1 in the SI). Our previous BOMD results\n5\nindicate the PW91-based ground state PES already described well the adiabatic VET\ndynamics [J. Phys. Chem. Lett. 10, 5969 (2019)] Altogether, these two commonly-used\nfunctionals representing two opposite limiting cases are sufficient for exploring how\nadiabatic and nonadiabatic energy transfer dynamics are influenced by the PESs (or\ndensity functional).\nA critical finding is that the accessibility to the barrier (state-crossing) region\nlargely controls both adiabatic and nonadiabatic VET dynamics. This means that the\nextent of both adiabatic and nonadiabatic energy transfer will vary consistently, and it\nis nearly infeasible to arbitrarily adjust their contributions to make a coincidental\nagreement with experiments under different conditions. As a result, the success of our\n‘golden combination’ method lies in its accurate treatment of both adiabatic and\nnonadiabatic energy transfer channels, rather than being ‘a fortuitous case’ or ‘an error\ncancellation’.\nTesting against more density functionals and identifying an ideal one is beyond\nthe scope of this work. However, considering the significant influence of the\ndissociation barrier on both adiabatic and nonadiabatic energy transfer, we calculated\nthe barrier height of NO on Au(111) using several other density functionals, including\na van der Waals corrected functional (vdW-DF), a meta-GGA functional (MN15-L) and\na meta-GGA functional plus van der Waals correction (MN15-L-rVV10), available in\nCP2K. The barrier heights predicted by these functionals are 2.77 eV, 3.06 eV, 2.86 eV,\nrespectively, which are close to that obtained with PW91. This is quite encouraging and\nindicating that the good performance of the PW91-based PESs is likely not accidental.\nThese additional barrier height tests have been added to Table S1 in the SI.\nChanges in revised manuscript: We changed the statement slightly to be more\nprecise on page 19 in the revised manuscript:\n“Given its significantly improved agreement with experimental results across\ndiverse conditions compared to previous models—unlikely due to significant error\ncancellations—the current PW91-based IESH model represents the state-of-the-art\ntheoretical characterization of molecule-surface scattering processes involving strong\nnonadiabatic effects.”\nComment 4: So, having found this ‘golden combination’ for NO on Au(111) is by\nitself important [indeed the research is also state-of-the-art] and deserves a publication\nin decent journal, it is however not a combination that includes a newly developed\n[general] method, neither is it a finding of new physics as such for this specific system\n6\nor a class of systems. I recommend therefore submitting this paper to another journal\nwhere then also the room of a more in-depth discussion is possible regarding the sources\nof errors that do remain.\nAuthor’s reply: We thank the reviewer for recognizing that our research is indeed\nstate-of-the-art. We would like to emphasize that the combination of CDFT, EANN,\nand IESH is certainly a general framework, which can be applied to study nonadiabatic\ndynamics in other gas-surface systems with electron transfer characteristics (as we will\nshow below for HCl/Au(111)). And IESH itself can be replaced by other nonadiabatic\ndynamical methods such as broadened classical master equation (BCME).\nAdditionally, several key experimental observations that have never been well\nunderstood before are now explained by our work, providing new physical insights into\nthis class of systems. For example, the variation in the steric effect of vibrational\nrelaxation is captured here by theory for the first time and is interpreted as a change in\norientational dependence of metal-to-molecule electron transfer behavior with the\nelongation of the molecular bond length. This new understanding would not be possible\nwithout an analysis of accurate diabatic PESs. Moreover, experimentalists have\npreviously speculated several possible mechanisms for the vibrational-to-translational\nenergy transfer of NO(v=3). Our calculations reveal that molecular vibration couples\ni\ndirectly to translation during the NO scattering, while the coupling between molecular\ntranslation and surface electrons is negligible. This identifies a direct vibrational-to-\ntranslational coupling, rather than an indirect electron-mediated coupling, as supposed\nin experiment. These findings would significantly advance our general understanding\nof energy transfer mechanisms at molecule-metal interfaces. As a result, we believe that\nthe general applicability of our method and new physical insights offered by our work\nsupport its publication to Nature Communications.\nChanges in revised manuscript: We further strengthened the discussion of the\nnew physical insights on page 23 in the revised manuscript:\n“Our simulations also provide new interpretations for some experimental\nobservations that were not well understood. For instance, the dependence of the steric\neffect on the initial vibrational state was attributed to the change of the orientational\ndependence of metal-to-molecule electron transfer behavior as the molecule stretches.\nAdditionally, our results indicate that during NO scattering, vibrational energy can\ncouple to both surface electrons and molecular translation, whereas the direct coupling\nbetween molecular translation and surface electrons is negligible.”\nGeneral comments: This remarkable article addresses a long standing and central\nissue in the dynamics of molecular interactions at surfaces – what is the relative\nimportance of effects that are within the Born Oppenheimer Approximation, i.e.,\nelectronically adiabatic effects and those that involve electronically nonadiabatic\ntransitions. They address this question by using state of the art theoretical methods to\ndo calculations modeling the rich trove of experimental work on NO vibrationally\ninelastic scattering from metal surfaces and compare models that that include or exclude\nnonadiabatic transitions but otherwise are on an equal footing with respect to the\npotential energy surface and dimensionality of the problem. The agreement of the best\nmodels to a broad range of experimental data as stunning and represents a tremendous\nadvance. Particularly remarkable is the ability of this work to reproduce the subtle\nresults on the initially molecular orientation dependance of vibrational energy transfer\nand its relation to crossing of the seam between different states on the potential energy\nsurface. The authors clearly demonstrate what they state in the title: both adiabatic and\nnonadiabatic effects matter.\nThese findings represent a dramatic advance in our understanding and ability to\ncalculate vibrationally energy transfer including nonadiabatic electronic effects that\nshould be of broad interest.\nI recommend publication in Nature Communications after the authors address the\nremarks given in the remainder of this review.\nAuthor’s reply: We sincerely thank Reviewer #2 for their encouraging feedback\nand thoughtful summary of the key contributions of our work. In response to the\nreviewer’s specific suggestions, we have carefully addressed all comments below and\nimproved the clarity and completeness of our discussion.\nComment 1: The abstract needs work.\nIt is OK up to the 4th sentence where the authors write “Our simulations reproduce,\nfor the first time, most experimental observations on vibrational relaxation and\nexcitation of NO molecules under varying initial conditions and clearly elaborate the\nrespective adiabatic and nonadiabatic contributions”.\nThis is followed by “This model identifies the direct adiabatic vibration-to-\ntranslation coupling as the dominant role of translation, while excluding significant\ntranslation-to-electron nonadiabatic coupling”.\nTaken in the context of the previous sentence the reader will likely think that that\n8\nthis sentence is providing elaboration on the respective adiabatic and nonadiabatic\ncontributions. Thus the sentence could be understood as saying that direct vibration-to-\ntranslation coupling is the dominant mechanism in vibrational relaxation. I think with\nthis sentence the authors are referring to the results on the transfer of vibrational energy\nto translation in vibrational deexcitation, which they show arises in an adiabatic model\nwith no participation of ehp coupling to translational motion of the NO molecule. But\nthis is not conveyed by saying “the model identifies the direct adiabatic vibration-to-\ntranslation coupling as the dominant role of translation.” Translation has many roles\nincluding pushing the NO to the seam in the PES, exciting phonons in the collision, etc.\nI suggest the authors add a sentence or two after “and clearly elaborate the\nrespective adiabatic and nonadiabatic contributions” to give some elaboration. For\nexample they might say something like vibrational relaxation of highly vibrationally\nexcited NO has comparable contributions from adiabatic and nonadiabatic effects. Then\nthey could say the coupling of the energy lost in vibrational deexcitation is\npredominantly adiabatic and is not significantly affected by nonadiabatic translation-\nto-electron coupling.\nAuthor’s reply: We appreciate the reviewer’s insightful and helpful suggestion.\nChanges in revised manuscript: To avoid the misleading expression, we have\nreformulated the abstract following the reviewer’s suggestion.\n“…… Our simulations reproduce, for the first time, most experimental\nobservations on vibrational relaxation and excitation of NO molecules under varying\ninitial conditions. We demonstrate comparable adiabatic and nonadiabatic contributions\nto the vibrational relaxation from initially low and high vibrational states. This model\nalso identifies that translational motion could couple to molecular vibration directly,\nwhile the translation-to-electron nonadiabatic coupling is not significant.”\nComment 2: Comparison of full-dimensional and two-dimensional dynamics.\nThis section would benefit from some introductory text to help the reader\nunderstand the motivations and questions that are being asked. There seem to be two\nmotivations: (1) to see how important multidimensional effects are and (2) to used 2D\ncalculations to validate IESH since a full quantum method is possible in 2D but not\ncurrently feasible in FD. Starting the section with some text to guide the reader rather\nthe starting with introducing a figure would help.\nAuthor’s reply: We thank the reviewer for this nice suggestion.\nChanges in revised manuscript: To clarify our motivation of comparing full-\n9\ndimensional and two-dimensional results, we added some introductory sentences on\npage 7 in the revised manuscript.\n“We first evaluate the accuracy of IESH by benchmarking IESH with a fully\nquantum method, namely hierarchical equations of motion (HEOM)40, 48. Since FD\nHEOM simulations are currently unfeasible, we compare IESH and HEOM results\nbased on a 2D model. As shown in Figs. 1a-c, IESH agrees well with HEOM, although\nthe latter predicts slightly stronger vibrational inelasticity than the former in the case of\nNO(v=3). …… To demonstrate the necessity of encompassing all DOFs in studying\ni\nthis energy transfer process, we also compare IESH results using 2D and FD\nmodels. ……”\nComment 3: It might be helpful to explain why was the revPBE functional chosen\nfor this work.\nAuthor’s reply: We chose the revPBE functional for two major reasons. First,\nIESH and HEOM results based on 2D PESs fitted from dozens of data points using the\nrevPBE functional were available, thus we can conveniently compare these results with\nour FD-IESH results to validate the accuracy of IESH and reveal the importance of\nmultidimensional effects. Second, the adsorption energy of NO on Au(111) predicted\nby revPBE is -0.14 eV, in reasonably good agreement with the experimental value (-\n0.24 eV, [Surf. Sci. 650, 11 (2016)]).\nChanges in revised manuscript: This part is added in the Methods section on\npage 25 of the revised manuscript.\n“We first chose the revPBE functional67 for two main reasons. First, IESH and\nHEOM results based on 2D (only includes the N-O bond length (r ) and the molecular\nNO\nheight on the surface (Z ), as defined in Fig. S1) PESs fitted from dozens of data points\nNO\nusing the revPBE functional were available40. This allowed us to conveniently compare\nthese results with our FD-IESH results to validate the accuracy of IESH and highlight\nthe importance of multidimensional effects. Second, the revPBE functional predicts a\nbinding energy of NO with Au(111) that agrees with the experiment value68 (see Table\nS1).”\nComment 4: Influence of PES on VET dynamics\nThis section could benefit from clearer descriptions. For example I can not\nunderstand what is meant by “while the MDEF(ODF) model further increases the\nprobability of single-quantum vibrational relaxation (i.e. ∆v=1)”. Please refer to which\n10\npanel is referred to.\nAuthor’s reply: We apologize for the ambiguous phrasing. As shown in Figs. 1d-\nf, the MDEF(ODF) model increases the probability of single-quantum vibrational\nrelaxation (i.e., ∆v=1) compared to BOMD on the PW91-based PES (see the\ncomparison between deep blue triangles and green rhombuses), resulting in better\nagreement with experimental observations. However, for multi-quantum vibrational\nrelaxation processes (i.e., ∆v > 1), the improvement by MDEF(ODF) is limited.\nChanges in revised manuscript: We have revised the relevant sentences on pages\n8-9 of the revised manuscript accordingly:\n“As discussed previously32, 33, 35, PW91-based BOMD results already predict\nremarkable vibrational relaxation. The MDEF(ODF) model further increases the\nprobability of single-quantum vibrational relaxation (i.e., ∆v = 1) compared to BOMD\non the PW91-based PESs (Figs. 1d-1f), leading to improved agreement with\nexperimental data. However, it falls short of improving the multi-quantum vibrational\nrelaxation (i.e., ∆v > 1)33.”\nComment 5: The statements about IESH get a little lost by being tacked on at the\nend of a paragraph focused on other things. The agreement with IESH should be\nemphasized more, for example by starting a new paragraph saying at the outset it gives\nthe best agreement with experiment to date.\nAuthor’s reply: We thank the reviewer for the suggestion of emphasizing the\nperformance of IESH on the PW91-based PESs.\nChanges in revised manuscript: Following the reviewer’s suggestion, we start a\nnew paragraph to describe IESH results on page 9 in the revised manuscript.\n“Encouragingly, the IESH model, which treats both electron transfer and EHP\nexcitations on equal footing, achieves the best agreement with experimental\nobservations across all initial conditions to date. Quantitatively, not only does the\nprobability of the two-quanta relaxation (v = 3 → v = 1) increase remarkably, but the\ni f\noverall vibrational state distributions of NO(v = 11 and 16) also shift to lower states,\ni\nboth of which align more closely with the experimental distributions. ……”\nGeneral comments: The manuscript addresses the study of a rather complex\nsystem, which remains the subject of ongoing debate: the reaction of NO with Au(111).\nIt employs a detailed methodology aimed at accurately reproducing the involved\nphysical phenomena, and the analysis is thorough and significant.\nAuthor’s reply: We thank the Reviewer #3 for the generally positive assessment\nof our work.\nComment 1: However, the title is unclear, as it implies that adiabatic and non-\nadiabatic channels can be treated independently, which seems incorrect. These channels\ncannot be fully separated. Calculations can be entirely adiabatic, but what does it mean\nto be completely non-adiabatic? Furthermore, while this work emphasizes the\nimportance of non-adiabatic effects, this issue was already raised in Tully and Co.’s\nseminal article [Science 326, 5954 (2009)].\nAuthor’s reply: We understand the referee’s concern. We note first that adiabatic\nand nonadiabatic channels cannot be distinguished in experiment at all. This is precisely\nwhy theory becomes essential to reveal the real energy transfer mechanism. We also\nagree that adiabatic and nonadiabatic channels cannot be completely separated in theory.\nBOMD calculations are entirely adiabatic, while in MDEF or IESH calculations, the\nadiabatic energy transfer channel remains open. Fortunately, we can compare results\nfrom BOMD and IESH simulations based on the same set of PESs to approximately\nestimate the relative contributions from both channels. For example, as shown in Figs.\n1d-f, our PW91-based BOMD results capture about half of the experimentally observed\nvibrational energy losses for three different initial vibrational states. The corresponding\nIESH results predict more substantial energy losses and are much closer to experimental\ndata. The difference between BOMD and IESH outcomes clearly demonstrates the\nsignificance of nonadiabatic effects, although a quantitative separation of the individual\ncontribution of the adiabatic and nonadiabatic channel is impossible, as the scattering\npath is altered when additional electronic degrees of freedom (DOFs) are involved. In\nour manuscript, we never mention anywhere that the energy transfer process is\ncompletely nonadiabatic.\nOur title only states that both adiabatic and nonadiabatic channels matter in\nmolecule-surface energy transfer processes, which is an important fact that has not been\nwell recognized in this community. All previous models (and experimental studies) to\nsome extent overemphasize the nonadiabatic character of such systems, where the\n12\nmetal-to-molecule electron transfer seems to be the only important thing. Our study,\nfor the first time, quantitatively reveals the significance of the adiabatic energy transfer\nchannel and the influence of the landscape of the barrier region on both adiabatic and\nnonadiabatic energy transfer. In this respect, our conclusion is not to emphasize the\nimportance of nonadiabatic effects alone. Instead, we find that both adiabatic and\nnonadiabatic channels play essential roles in vibrational energy transfer dynamics. This\nfinding is completely different from the issue addressed by Tully and coworkers’\nseminal article [Science 326, 5954 (2009)].\nWe hope this clarification validate the use of the current title.\nChanges in revised manuscript: Inspired by the reviewer’s comment, we added\nsome discussions to compare our results with Tully’s work on page 9 in the revised\nmanuscript:\n“…… By comparing final vibrational state distributions predicted by BOMD and\nIESH, we conclude that adiabatic and nonadiabatic contributions to the vibrational\nrelaxation with low and high initial vibrational states are comparable. In contrast, in\nTully’s work21, they stated that the nonadiabatic channel dominates the VET for the\nNO/Au(111) scattering, which differs a lot from our conclusion, mainly resulted from\ntheir less accurate ground and charge-transfer state PESs.”\nComment 2: The comparison between the new theoretical results and\nexperimental data is outstanding and underscores the importance of performing\ndynamics on a highly accurate potential energy surface. Nonetheless, there is some\ndebate regarding the statement that the agreement cannot be due to error compensation,\nespecially considering that the potential energy surface was calculated using the PW91\nfunctional.\nAuthor’s reply: We sincerely thank the reviewer for approving the outstanding\nagreement between our theoretical results and experimental data. We argue that the\ngood agreement of our model with experiments under various conditions is not likely\ndue to significant error cancellations for several reasons. First, as discussed in our\nmanuscript (pages 9-11), the accessibility to the barrier (or state-crossing) region\npredominantly controls both adiabatic and nonadiabatic VET dynamics. Consequently,\nthe extent of adiabatic and nonadiabatic energy transfer will vary in a consistent manner,\nmaking it hardly feasible to arbitrarily adjust their contributions to achieve coincidental\nagreement with experiments under different conditions. Moreover, our additional tests\non the dissociation barrier height of NO on Au(111) also suggested that the barrier\n13\nheight obtained with PW91 is likely reproduced by other more sophisticated functionals,\nsupporting that its good performance is not a fortuitous coincidence (See Table S1 and\nour response to Reviewer #1, Comment 3). Second, the IESH method shows good\nagreement with a fully quantum method HEOM (see Figs. 1a-c), validating the\naccuracy of IESH and excluding significant errors from the dynamical method. As a\nresult, we tend to conclude that the success of our method mainly lies in its accurate\ntreatment of both adiabatic and nonadiabatic energy transfer channels, rather than being\na result of ‘error cancellation’.\nChanges in revised manuscript: We changed the statement slightly to be more\nprecise on page 19 in the revised manuscript:\n“Given its significantly improved agreement with experimental results across\ndiverse conditions compared to previous models—unlikely due to significant error\ncancellations—the current PW91-based IESH model represents the state-of-the-art\ntheoretical characterization of molecule-surface scattering processes involving strong\nnonadiabatic effects.”\nComment 3: Overall, the results are compelling and warrant publication. The\nmethodology includes all necessary components to reproduce the phenomena, and the\nfindings are meaningful. However, the conclusions may overstate the significance of\nthe work. It might be more appropriate for publication in a specialized journal.\nOn the other hand, the impact of this research could be enhanced by demonstrating\nthat the same methodology can be applied to other complex systems, such as\nHCl/(Au(111)). This would showcase the broader applicability and relevance of the\napproach.\nAuthor’s reply: We sincerely appreciate your helpful suggestion. We have\nreemphasized the significance of our work in our response to your Comment 1. To\nfurther validate the broad applicability of the proposed CDFT+EANN+IESH scheme,\nwe have extended this methodology to HCl/Au(111) following the reviewer’s\nsuggestion, and present some preliminary results here for review only. To this end, we\nhave first trained six-dimensional neural network PESs for the ground, neutral and\nanionic states of HCl/Au(111), including all molecular DOFs, based on several hundred\n(constrained) DFT data points using the revPBE functional (previous studies indicated\nthat this system is highly activated and requires a strongly repulsive density functional).\nFig. R2 shows the charge-transfer states and the ground state of this system. As\nexpected, electron transfer from Au(111) to HCl is much more difficult than to NO,\n14\nbecause of the more negative electron affinity of the HCl molecule.\nA key experimental observation for HCl/Au(111) is a transition of the vibrational\nexcitation mechanism from adiabatic to nonadiabatic as the surface temperature (T )\ns\nincreases for HCl(v = 0) scattering from the Au(111) surface [Phys. Rev. Lett. 98,\ni\n237601 (2007); J. Chem. Phys. 145, 054709 (2016)]. At low T , the vibrational\ns\nexcitation of HCl(v = 0) is likely dominated by the mechanical coupling between\ni\nmolecular vibration and translation, while at high T , thermally excited electron-hole-\ns\npairs will dominate the vibrational excitation. Encouragingly, our preliminary IESH\nresults with electronic temperatures ranging from 273 K to 1073 K at E=0.86 eV,\ni\nreproduced this experimental trend well, as shown in Fig. R3, although the absolute\nvibrational excitation probabilities are overestimated due possibly to the neglect of the\nlattice motion. Anyhow, these results validate the generalizability of the methodology\nand further strengthen the significance of our work.\nFig. R2. One-dimensional cuts of the PESs for the HCl/Au(111) system. (a) Energy\ncurves of the ground (E ), neutral (U ) and anionic (U ) states for the HCl/Au(111)\ng 0 1\nsystem as a function of the HCl molecular height above the surface (Z ), with the\nHCl\nbond length of HCl (r ) fixed at 1.28 Å (equilibrium bond length of the HCl molecule).\nHCl\nThe HCl molecule is oriented perpendicularly to the hcp site of the Au(111) surface,\nwith the H atom pointing toward the surface. (b) Similar to (a), but with r = 1.52 Å\nHCl\n(approximately the classical outer turning point of HCl(v = 1)). (c) Energy curves as a\nfunction of r with Z fixed at 2.50 Å. The energy zero is referenced to E of a free\nHCl HCl g\nHCl molecule far from the surface.\npredicted vibrational excitation (v = 0→1) probabilities during HCl(v = 0, E = 0.86 eV)\ni i\nscattering from Au(111) at various electronic temperatures with experimental results [J.\nChem. Phys. 145, 054709 (2016)].\n16"}, {"round": 2, "reviewer_comments": {"reviewer_1": "1. The authors have tried to address my point regarding another source of error in their work. They rightfully now include the\nPES and the used NN as a source, however, in determining, or elaborating on its accuracy, ambiguity is being introduced in\ntheir wording: \"To further validate the accuracy of the trained PESs, we have performed additional tests by randomly\nexcluding 10% of data points (~200 points) from the original data set and retraining the NN PESs.\".\nThis could mean two things. a) the original data set was [randomly] reduced by 10% and then a NN PES was trained from\nscratch on the reduced data set and the now reported RSMEs are based on the data now NOT included into the training, or\nb) the fully trained NN that was trained on the complete data set is used as a starting point and 10% of the original data set is\n[randomly] cut out and the already trained NN is then given a few more training epochs and when it then again converges,\nthe RSMEs are computed for the part of the data set that was cut out.\nA proper test would be doing a) where the NN is tested on data that has NOT anywhere or EVER been used during training.\nIt would be no surprise to me that when b) was done, there are hardly any changes as the grunt of the training of the NN was\ndone whilst using the whole data set.\nThe wording needs to be clear here if a) was done. If b) was done, the test is just no reflection of the quality of the NN at all\n[only its converging is tested, not to WHAT it actually converges too when started from scratch].\n2. The authors give a lengthy reply to confirm that the use of any DFT functional is a source of error in their work and point\nout that this is a typical situation in the field. Even for state-of-the-art studies. I am not disagreeing with that. It also precisely\nunderlines my original statement about the found 'golden combination' to be probably a fortuitous case and far from being a\ngeneral find or a \"breakthrough\".\nLet me try and explain it in an car analogy. Suppose Ferrari has overhauled their F1 car and sourced many parts from many\nvendors and factories, and then has a very fast lap on one particular track only. However, when one part is exchanged or a\nrace is run on another track, their performances drop significantly and they again end up among the rest of the competition.\nNo one in his right mind would then make the claim that Ferrari has made a genuine breakthrough in F1 racing!\nYet this is what the authors still seem to do.\nIn their rebuttal a lengthy section regarding the difficulties of DFT is given in the field [all examples given, I am also aware\nof]. However no new statement, fact or reason is given why their RESULTS are to be considered a breakthrough. Nor is\nthere a new fact or argument given on why this inherent dependency of their RESULTS on having everything carefully\naligned into the 'golden combination' of elements [one of these element is thus the specific DFT functional being required],\nis to be considered not to be a weak point. They simply only state it should be considered to be a strong point after a lengthy\npiece of text on how difficult things are with DFT functionals in the field?\nI therefore personally remain unconvinced that this work deserves a publication in Nature Communication.\nI still do think the work itself would be of interest to the community in a more suitable journal that also allows the proper\ntreatment [and not an ad-hoc fix perhaps aimed to only please a referee or editor?] of my worries regarding the NN PES in\nthe point 1 above.", "reviewer_2": "I thank the authors for their thoughtful consideration of the points made in my review and the other review. The manuscript\nhas been improved by their revisions in response to these suggestions.\nI recommend publication of manuscript NCOMS-25-2791A in its revised form in Nature Communications.\nDaniel Auerbach", "reviewer_3": "In my opinion, the authors have effectively addressed all of my concerns, along with those raised by the other reviewers. The\nrevisions made to the manuscript have enhanced the clarity and comprehensibility of its main message, emphasizing its\nsignificance. Consequently, I recommend the publication of this manuscript in its current form."}, "author_response": "Response to Reviewers\nNCOMMS-25-27914A: First-principles understanding of vibrational energy\ntransfer in molecule-surface scattering: both adiabatic and nonadiabatic\nchannels matter\nReviewer Comments: Blue\nAuthor Responses: Black\nRevisions: Red\nComments from Reviewer #1 and our detailed responses\nComment 1: The authors have tried to address my point regarding another source\nof error in their work. They rightfully now include the PES and the used NN as a source,\nhowever, in determining, or elaborating on its accuracy, ambiguity is being introduced\nin their wording: \"To further validate the accuracy of the trained PESs, we have\nperformed additional tests by randomly excluding 10% of data points (~200 points)\nfrom the original data set and retraining the NN PESs.\".\nThis could mean two things. a) the original data set was [randomly] reduced by\n10% and then a NN PES was trained from scratch on the reduced data set and the now\nreported RSMEs are based on the data now NOT included into the training, or b) the\nfully trained NN that was trained on the complete data set is used as a starting point and\n10% of the original data set is [randomly] cut out and the already trained NN is then\ngiven a few more training epochs and when it then again converges, the RSMEs are\ncomputed for the part of the data set that was cut out.\nA proper test would be doing a) where the NN is tested on data that has NOT\nanywhere or EVER been used during training. It would be no surprise to me that when\nb) was done, there are hardly any changes as the grunt of the training of the NN was\ndone whilst using the whole data set.\nThe wording needs to be clear here if a) was done. If b) was done, the test is just\nno reflection of the quality of the NN at all [only its converging is tested, not to WHAT\nit actually converges too when started from scratch].\nAuthor’s reply: We thank the reviewer for pointing out the possible ambiguity.\nWe clarify that workflow (a), as described by the reviewer, was done to validate the\naccuracy of the trained PESs. Specifically, we randomly excluded 10% of the original\ndataset (~200 data points) and then trained new PESs from scratch using only the\nremaining 90% of the data points. These newly trained PESs were then used to compute\n1\nthe test RMSEs on the 10% of data points that were completely excluded from the\ntraining process, and perform nonadiabatic simulations under the condition of NO(v =\ni\n16, E = 0.52 eV). We believe that these tests are sufficient to validate the accuracy of\ni\nthe NN PESs used in our work.\nChanges in revised Supplementary Information (SI): To clarify our approach,\nwe have revised the description on page S14 of the SI as follows:\n“To validate the convergence of the used PESs, we performed additional tests. We\nrandomly excluded 10% of the data points from the initial dataset and trained new PESs\nfrom scratch using only the remaining 90% of the data points. The test RMSEs\ncomputed for the excluded 10% of data points using these newly trained PESs were\n0.15, 0.56, and 0.57 meV/atom for the ground, neutral, and anionic PESs,\nrespectively—values comparable to the total RMSEs of the original PESs. ……”\nComment 2: The authors give a lengthy reply to confirm that the use of any DFT\nfunctional is a source of error in their work and point out that this is a typical situation\nin the field. Even for state-of-the-art studies. I am not disagreeing with that. It also\nprecisely underlines my original statement about the found 'golden combination' to be\nprobably a fortuitous case and far from being a general find or a \"breakthrough\".\nLet me try and explain it in an car analogy. Suppose Ferrari has overhauled their\nF1 car and sourced many parts from many vendors and factories, and then has a very\nfast lap on one particular track only. However, when one part is exchanged or a race is\nrun on another track, their performances drop significantly and they again end up\namong the rest of the competition.\nNo one in his right mind would then make the claim that Ferrari has made a\ngenuine breakthrough in F1 racing!\nYet this is what the authors still seem to do.\nIn their rebuttal a lengthy section regarding the difficulties of DFT is given in the\nfield [all examples given, I am also aware of]. However no new statement, fact or reason\nis given why their RESULTS are to be considered a breakthrough. Nor is there a new\nfact or argument given on why this inherent dependency of their RESULTS on having\neverything carefully aligned into the 'golden combination' of elements [one of these\nelement is thus the specific DFT functional being required], is to be considered not to\nbe a weak point. They simply only state it should be considered to be a strong point\nafter a lengthy piece of text on how difficult things are with DFT functionals in the field?\nI therefore personally remain unconvinced that this work deserves a publication in\nI still do think the work itself would be of interest to the community in a more\nsuitable journal that also allows the proper treatment [and not an ad-hoc fix perhaps\naimed to only please a referee or editor?] of my worries regarding the NN PES in the\npoint 1 above.\nAuthor’s reply: Our previous rebuttal was necessarily detailed to\ncomprehensively address the referee’s Comments #2, #3, and #4, all of which are\nrelevant to the so-called ‘golden combination’. However, ONLY in one paragraph,\nwhich is not too lengthy, we listed a few examples to illustrate the intrinsic dependence\nof dynamical results on DFT functionals and, more generally, the fact that any\ntheoretical study of chemical dynamics relies on a combination of methods (e.g.,\nelectronic structure calculations, PES fitting, and nuclear or electron-nuclear coupled\ndynamics).\nTo make full-dimensional nonadiabatic dynamics simulations feasible for this\nchallenging system, our model integrates currently most adequate approaches,\nachieving an unprecedented level of accuracy when using the PW91 functional. We\nhave justified our method choices in our previous rebuttal in several ways. First, we\ngave a new fact that the barrier height obtained by PW91 (GGA) is close to those\npredicted by more sophisticated functionals (meta-GGA and vdW correction). Second,\nwe noted that the IESH method and a fully quantum HEOM method yield similar results\nin 2D cases. We further argued that because both adiabatic and nonadiabatic energy\ntransfer are controlled by the barrier (state-crossing) region, it is difficult to arbitrarily\nadjust their relative contributions to reproduce experiments under different conditions.\nAs a result, the general agreement between our results and diverse experimental data is\nmore likely due to our accurate treatment of both adiabatic and nonadiabatic energy\ntransfer channels, rather than a fortuitous combination of a specific DFT functional and\ndynamical theory.\nIn our view, on the premise of the good agreement between theory and experiment,\nwhat is more important for a theoretical work lies in elucidating detailed mechanisms\nbehind experimental observations that have not been well understood. In our previous\nrebuttal, therefore, we have devoted more effort to articulating why our work represents\na breakthrough in this field, meriting publication in Nature Communications rather than\nin a more specialized journal. Emphasizing again, through a first-principles, full-\ndimensional characterization of NO scattering from Au(111), our work provides\nconvincing evidence for the critical roles of both adiabatic and nonadiabatic VET\n3\npathways in this prototypical molecule-metal system, while previous studies mostly\noverlooked the adiabatic contribution. Moreover, our analysis clearly demonstrates how\nthe potential energy landscape influences both adiabatic and nonadiabatic VET\nprocesses (Figs. 1 and 2), which would not be possible without accurate CDFT-based,\nEANN-fitted diabatic PESs. These findings provide new fundamental insights into\nenergy transfer dynamics during molecular collisions on metal surfaces. Additionally,\nour results account for several experimental observations that had never been\nreproduced by theory before. For example, the dependence of the steric effect on the\ninitial vibrational state is attributed to the change in the orientational dependence of\nmetal-to-molecule electron transfer behavior as the molecule stretches. We further\nidentify that molecular vibration couples strongly to both surface electrons and\nmolecular translation, whereas direct coupling between molecular translation and\nsurface electrons is negligible (Figs. 4 and 5). These findings establish a complete\nmechanistic picture of energy transfer dynamics in this representative system and reveal\nthe interplay between competing channels at the molecule-metal interface. Overall, our\nwork not only reproduces a wide range of experimental data, but also provides valuable\ninsights into energy transfer dynamics in this challenging system and, more generally,\nat molecule-metal interfaces.\nLet us contextualize our breakthrough using the referee’s F1 analogy. From one\nperspective, each track could correspond to a certain experimental condition. The\nearlier model combining IESH with empirical RST PESs [Science 326, 829 (2009)]\nthen resembles a race car that has a very fast lap on a particular track only, but suffers\na dramatic performance loss on other tracks. Specifically, that model coincidentally\nreproduced experimental results for NO(v=15) scattering at low incidence energies,\ni\nwhile failed under many other experimental conditions [J. Chem. Phys. 140, 044701\n(2014), J. Phys. Chem. C 119, 3268 (2015)]. In contrast, our model achieves consistent\nagreement with most experimental data, and we have explained why this is the case and\ndiscussed the direction for further improvement (see Discussion section). In this\nanalogy, our model is the ONLY race car equipped with all essential components, fast\nenough on most tracks, and with clear directions for further improvement. This, in our\nview, represents a genuine breakthrough in the field.\nIn an alternative perspective, each track may correspond to a specific system, so\nthe referee is questioning the generalizability of our model to other systems. In this\nanalogy, the NO/Au(111) system represents one of the most challenging tracks, and our\nmodel resembles currently the ONLY car that successfully completes the race, setting\n4\na record time. Moreover, our results for CO/Au(111) [Phys. Rev. Lett. 133, 036203\n(2024)] and preliminary results for HCl/Au(111) (presented in our previous rebuttal),\neach based on a suitable density functional, demonstrated the general applicability of\nour proposed CDFT+EANN+IESH scheme. To the best of our knowledge, no\ncompeting schemes have achieved comparable performance across such diverse\nsystems. This indicates that our car can adapt to other tracks by adjusting only one\ncomponent upon tests (e.g., selecting an appropriate density functional for a given\nsystem), much like F1 teams optimize their cars for different tracks, while at present,\nno other competitors are even able to finish the race on these tracks.\nFrom either perspective, we think the referee’s analogy is not weakening but\nhelping show the excellence of our work. We hope our arguments above will convince\nthe referee that the combination of methods should not be detracted as a weak point,\nparticularly directed against to our work. In fact, ‘golden combination’ itself sounds\nlike that the combination is intentionally chosen to match experiments, which is an\ninappropriate assessment yet not the truth.\nComments from Reviewer #2 and our detailed responses\nGeneral comments: I thank the authors for their thoughtful consideration of the\npoints made in my review and the other review. The manuscript has been improved by\ntheir revisions in response to these suggestions.\nI recommend publication of manuscript NCOMS-25-2791A in its revised form in\nNature Communications.\nAuthor’s reply: We appreciate the referee’s recommendation.\nComments from Reviewer #3 and our detailed responses\nGeneral comments: In my opinion, the authors have effectively addressed all of\nmy concerns, along with those raised by the other reviewers. The revisions made to the\nmanuscript have enhanced the clarity and comprehensibility of its main message,\nemphasizing its significance. Consequently, I recommend the publication of this\nmanuscript in its current form.\nAuthor’s reply: We appreciate the referee’s recommendation.\n5"}]

Decision letter