An experimentally validated neural-network potential energy surface for H atoms on free-standing graphene in full dimensionality

We present a first principles-quality potential energy surface (PES) describing the inter-atomic forces for hydrogen atoms interacting with free-standing graphene. The PES is a high-dimensional neural network potential that has been parameterized to 75945 data points computed with density-functional theory employing the PBE-D2 functional. Improving over a previously published PES (Jiang et al., Science, 2019, 364, 379), this neural network exhibits a realistic physisorption well and achieves a 10-fold reduction in the RMS fitting error, which is 0.6 meV/atom. We used this PES to calculate about 1.5 million classical trajectories with carefully selected initial conditions to allow for direct comparison to results of H- and D-atom scattering experiments performed at incidence translational energy of 1.9 eV and a surface temperature of 300 K. The theoretically predicted scattering angular and energy loss distributions are in good agreement with experiment, despite the fact that the experiments employed graphene grown on Pt(111). The remaining discrepancies between experiment and theory are likely due to the influence of the Pt substrate only present in the experiment.


Introduction
H-atom chemisorption to graphene is relevant to hydrogen storage 1 , the catalytic formation of molecular hydrogen in the interstellar medium 2 and -because hydrogenation of graphene can induce a band-gap -two-dimensional semiconductor materials 3 . Recently, a full-dimensional PES was reported using first principles energies obtained from Embedded Mean-field Theory (EMFT) [4][5][6] to parameterize a second generation Reactive Empirical Bond Order (REBO) function 7 . Using classical and semiclassical dynamics calculations, qualitative agreement was obtained with H-atom scattering experiments carried out at incidence translational energies E i of 1.9 and 1 eV. Furthermore, the trajectories provided an atomic scale movie at the femtosecond time scale showing the formation of a covalent chemical bond 8 . The sticking probability could also be calculated using the REBO-EMFT PES and compared well with experiment at E i = 1 eV. This suggests that the REBO-EMFT PES is the best available representation of interatomic forces in the H/graphene system.
Despite the progress made in that work, two problems remained. First, the EMFT data was derived from a model of free-standing graphene, while the experiment was carried out on graphene that had been grown on Pt(111) 8 . To account for the influence of Pt in the simulations, a Lennard-Jones (LJ) interaction model with the Pt substrate was included for each atom in the graphene layer. This improved agreement with experiment, suggesting the influence of the substrate may be important. Unfortunately, it is unclear how to reparameterize the analytical REBO-EMFT PES from first-principles energies that include the Pt substrate. Therefore, the role of the substrate remains uncertain. The second problem concerns the fitting error (7 meV/atom) as the REBO function is not flexible enough to closely reproduce electronic structure data 8 . This complicates the evaluation of the quality of different electronic structure methods, since the fitting error can easily be larger than the energy differences between the methods being compared. Clearly, a full-dimensional firstprinciples PES where fitting errors are small and where the role of the Pt substrate is included would be a significantly better approach to this problem. For both of these problems a solution is offered by atomistic potentials employing machine learning (ML) methods.
In recent years, ML potentials have become a promising new approach to construct PESs of first-principles quality 9,10 . They have a uniquely flexible functional form that allows the accurate reproduction of reference data sets obtained in electronic structure calculations, without sacrificing the efficiency needed when they are repetitively evaluated in large-scale molecular dynamics (MD) simulations. ML potentials have been developed for many systems. These include free-standing 11 and multi-layer graphene 12 as well as graphite [12][13][14] and amorphous carbon 15 , which are closely related to this work. A frequently used type of machine learning potential suitable for large condensed systems is the high-dimensional neural network potential (HDNN-PES) method proposed by Behler and Parrinello in 2007 16 .
In this paper, we present the first HDNN-PES for H atoms interacting with free standing graphene, which we validate against data obtained from H and D scattering experiments using graphene grown on Pt(111), experiments that are similar to those recently reported elsewhere. 17 Compared to the REBO-EMFT PES 8 , we achieve substantially reduced fitting errors without sacrificing computational performance. Using molecular dynamics, we show that experimentally obtained H/D-atom energy loss and angular distributions are faithfully reproduced. We demonstrate the improvement represented by the HDNN-PES by comparing the new results to MD simulations done with the previously reported REBO-EMFT PES 8 . The remaining deviations between experiment and theory likely reflect the absence of the Pt substrate in our simulations; however, the influence of the substrate on the scattering distributions appears to be relatively small.

Experimental Methods
The experimental apparatus has been described in detail in Ref. 18. H/D-atoms are generated by photodissociating a supersonic molecular beam of hydrogen/deuterium iodide with a KrF excimer laser producing atoms with incidence energy of 1.92 eV. A small fraction of these atoms passes through two differential pumping stages, enter the ultra-high vacuum chamber and collide with the graphene sample grown in situ on a Pt(111) substrate. The sample is held on a six-axis manipulator, allowing variation of the incidence angle θ i . Recoiling atoms are excited to a long lived Rydberg state (n = 34) by two laser pulses at 121.5 nm and 365 nm via a two-step excitation. These neutral atoms travel 25 cm in a field-free region and pass a detector aperture before they are field-ionized and detected by a multi-channel plate detector. The arrival time is recorded by a multi-channel scalar. The rotatable detector allows data to be recorded at various scattering angles θ s . The graphene sample is epitaxially grown on a clean Pt(111) substrate by dosing ethylene (partial pressure 3 × 10 −8 mbar) at 700 • C for 15 minutes.

HDNN-PES
High-dimensional neural network potentials (HDNN-PESs) 16 have been the first type of ML potential enabling the simulations of large condensed systems. In this approach, the total potential energy E tot of the system is constructed as a sum of atomic energy contributions, depending on the local chemical environment defined by a cutoff radius R c , typically in the range between 6 and 10 Å. The positions of all neighboring atoms inside the cutoff spheres are described by sets of atom-centered many-body symmetry functions 19 . The resulting vector of symmetry function values for each atom represents a structural fingerprint that is used as input for an atomic neural network yielding atomic energy contribution E µ into the total energy (1). The functional forms of the symmetry functions ensure the necessary invariance of PES with respect to translations and rotations of the system as well as permutations of like atoms. The atomic neural networks are feed-forward neural networks and contain a large number of weight parameters, which serve as fitting parameters for the HDNN-PES. Each element in the system is modeled by a separate atomic neural network with a specific architecture and values of the weight parameters cal-culated once for each atom of the respective element in the system. The values of these parameters are determined in an iterative optimization process by minimizing the errors of the energies and forces for a reference data set of representative structures obtained from electronic structure calculations, typically densityfunctional theory. Additional structures that may be required in the reference set in regions of the PES that are not well sampled can be suggested by an automatic procedure employing a committee of HDNN-PESs and a comparison of predicted energies and forces 20 leading to a self-consistent and unbiased generation of the data set. Once a set of weight parameters accurately reproducing the reference data has been found, the PES undergoes a series of careful validation steps 21 . Then, the HDNN-PES is ready for applications. For all the details about the HDNN-PES method, the determination of the weight parameters and its validation procedures, the interested reader is referred to several recent reviews 21-23

Density Functional Theory Calculations
The Vienna ab initio simulation package (VASP) version 5.3.5 [24][25][26][27] has been employed for the reference electronic structure calculations to generate the training set for the HDNN-PES. Density functional theory (DFT) at the generalized gradient approximation (GGA) level of theory using the Perdew-Burke-Ernzerhof (PBE) 28 exchange-correlation functional with a planewave basis has been used in combination with Grimme D2 van der Waals (vdW) corrections 29 . We made use of the Projector Augmented Wave (PAW) 27,30 approach to model the core and valence electron interactions. The kinetic energy cutoff has been set to 400 eV. The Monkhorst-Pack scheme 31 with a 8 × 8 × 1 Γcentered k-point mesh for the 3 × 4 surface cell has been used to sample the surface Brillouin zone. With two atoms per primitive unit cell, the slab consists of 24 carbon atoms in total and is 8.55 Å × 7.40 Å in size (see Fig. 1). 3D periodic boundary conditions have been applied with 13 Åvacuum perpendicular to the graphene sheet to ensure that the periodic images of the surfaces are non-interacting and that hydrogen atoms can be included at a maximum separation of 6 Å from the surface. We included spin polarization in the electronic structure calculations, and partial occupations have been treated by applying the tetrahedron method with Blöchl corrections 30,32 using the default value of 0.2 eV as the smearing parameter. The threshold for the change in energy between iteration steps when relaxing the electronic degrees of freedom has been 10 −5 eV.

Generating the Reference Structures
The iterative procedure described in detail elsewhere 34 was used to generate the reference data. Briefly, step by step new DFT energies and forces are added for geometries where the HDNN-PES fit does not show the desired accuracy or covers the full configurational space. These geometries are identified by comparing the results of several generated HDNN-PESs with differing network structures. The reference data set is then extended with the additional data until convergence is reached. The initial data set consists of energies and forces obtained 3.4 Construction of the Neural Network Potential from DFT calculations for about 6 × 10 4 reference configurations, which were picked up from: (i) ab initio molecular dynamics (AIMD) trajectories simulating H-atom scattering from graphene at incidence energy of 1.9 eV and incidence angles of 34 • and 52 • at surface temperatures of 300 K and 600 K; (ii) geometries close to the minimum energy path to adsorption, where the Hatom was put at the lateral position of the C-atom and the zcoordinates were varied over a range of −0.8 Å ≤ H z ≤ 5.8 Å and −0.8 Å ≤ C z ≤ 1.0 Å, respectively, with 0.025 Å step and without structures with r CH < 0.6 Å, whereas the remaining C-atoms were kept at their equilibrium positions; (iii) graphene geometries chosen randomly from an AIMD trajectory thermalized at 300 K with H-atom over the surface. The position of the H-atom is chosen randomly over the whole simulation cell where the z-coordinate ranges from 1 to 6 Å. The configuration with a C-H distance of 6 Å and a fully relaxed graphene surface was used as the asymptotic energy reference. This structure is our energy zero point. The HDNN-PESs fitted to the initial reference data set were then improved on the set of about 1.5 × 10 4 configurations obtained from MD simulations of H-atom scattering from a graphene sheet at incidence energy of 1.9 eV in the wide range of incidence angles (from 0 • to 90 • in 10 • step) as well as at incidence energy of 6 eV and normal incidence angle with surface temperatures of 0 K and 600 K starting over high-symmetry sites shown in Fig. 1. Moreover, we also trained the HDNN-PESs on the configurations taken from equilibrium MD simulations of the graphene surface in the wide range of temperatures from 0 K to 2000 K. The high-temperature configurations are useful, since the surface can be heated locally in the neighborhood of the collision site.
In total, the final HDNN-PES was trained on the reference data set of 75 945 configurations.

Construction of the Neural Network Potential
The HDNN-PES has been constructed using the RuNNer 21,22,35 code. The atomic neural network's architecture consists of two hidden layers with 15 neurons per layer providing the energies both for hydrogen and carbon atoms. The parameters of symmetry functions 19 are listed in Table 1 of SI. The symmetry function values have been rescaled to the range from 0 to 1. Randomly selected 90% of the reference data were used to train the NN, whereas the remaining 10% were used as an independent test set to validate the fit and to check for overfitting. The NN weight parameters were determined from the DFT energies and forces employing the adaptive global extended Kalman filter 36 . The initial values of the weight parameters have been chosen randomly in the interval from −1 to 1. For the weights, a preconditioning scheme was applied to reduce the initial root-mean-square error (RMSE) 21 . The training data in each of the 200 iterations (epochs) of the fit were presented in a random order to reduce the probability of getting trapped in local minima.

Molecular Dynamics Simulations Details
MD simulations of H-atom scattering from graphene were performed using MDT2 code 37 developed to study atomic scattering from various surfaces 8,38,39 . The RuNNer subroutines implementing the HDNN-PES providing the energies and forces were integrated into the MDT2 code. All the results shown in this paper have been obtained from MD simulations carried out using the RuNNer-MDT2 interface.
MD simulations of the H/D scattering from graphene have been carried out in the NVE ensemble using the standard velocity Verlet algorithm 40,41 with a time step of 0.1 fs. The trajectories were started with an H-atom randomly put at height of 3.5 Å over the surface and were terminated either when the scattered atom distance from the surface became larger than 3.6 Å or when the trajectory duration exceeded 200 fs. The initial geometry for the graphene layer was randomly selected from 1000 configurations obtained after the equilibration of the surface at 300 K with Andersen thermostat 42,43 . Those configurations were extracted from a single 100 ps-long trajectory with a period of 100 fs. In the experiment, graphene is not a single crystal but a composition of two equally abundant orientational domains. The two domains have a rotational distribution with a Gaussian width of 5 • and they are rotated by 27 • with respect to each other 8 . This results in a H-atom velocity vector that is oriented symmetrically with respect to the two domains. To achieve scattering conditions comparable to experiment, the simulations have been carried out by averaging over two domains with incidence azimuth φ i = ±13.5 • , where zero for azimuth angle is aligned with a C=C bond in graphene. In total, ∼1.5 million trajectories have been carried out for different incidence angles and isotopes (hydrogen and deuterium). The exact numbers of trajectories for the different conditions can be found in Table 2 of the SI. Fig. 2 shows a two dimensional cut through the converged HDNN-PES developed in this work reflecting structures near the minimum energy path to chemisorption, where the H-atom approaches directly above a C-atom. A physisorption well can be seen at large H z and a deeper chemisorption well at small H z with C z ≈ 0.4 Å. The minimum energy path to chemisorption in- volves both degrees of freedom, demonstrating that the C-atom is partially re-hybridized from sp 2 to sp 3 at the transition state. The depth of the physisorption well of the HDNN-PES has a depth of 9 meV. This compares well with the DFT energy calculated at this configuration 22 meV. The global physisorption minimum has the H-atom centered over the 6-ring. Here, the HDNN-PES gives the well depth of 11 meV at an H z = 2.7 Å. This compares reasonably well with the experimentally determined physisorption well depth (40 meV) 44 and a correlated, counterpoise corrected wave function calculation of the hydrogen-coronene system, which found the minimum at H z = 2.93 Å 45 . The previous REBO-EMFT PES had no physisorption well.

Results
The chemisorption well depth of the HDNN-PES (657 meV) also compares well with DFT (676 meV) but is deeper than that of the REBO-EMFT PES (610 meV). Furthermore, the DFT barrier (160 meV) is reproduced well by the HDNN-PES (172 meV) but is lower than that of the REBO-EMFT PES (260 meV).
The improved quality of the HDNN-PES in comparison to the REBO-EMFT PES is due both to the use of a dispersion corrected functional as well as to reduced fitting error. The RMSE fitting error of the REBO function to the DFT training data was reported to be ≈ 7 meV/atom 8 ; furthermore, the REBO function cannot represent a physisorption well. The flexibility of the neural network-the RMSE for the HDNN-PES is ≈ 0.6 meV/atom for energies in training and test set and ≈ 90 meV/Å for forces in training and test set, respectively-easily leads to a physically realistic physisorption well and a more accurate representation of the DFT energies and forces. Fig. 3 shows the fitting error to the DFT energies graphically. While the errors are not randomly distributed, there is no reason to suspect systematic problems with the PES over the energy range of 10 eV.  AIMD. The trajectories correspond to the same initial conditions and are typical of those that will be compared to experiment below. Fig. 4 shows the potential energy change along the trajectory while Fig. 5 shows the H atom motion in the trajectory in a perspective drawing. The trajectories obtained with these two approaches are nearly identical-note that if the fitting were perfect they would be identical. We now turn to the question: How well is the experiment reproduced by classical MD on the new HDNN-PES.
Figs. 6 and 7 show comparisons between experiment and theory for H and D scattering from graphene, respectively. In both Fig. 5 The same two trajectories as in Fig. 4-AIMD (red) and HDNN-PES (blue). The H-atom's initial position is shown as a cyan colored ball. The divergence between the two trajectories is due to residual error in the NN fit to the DFT data. A "side view" of the trajectories can be found in the SI.
figures, panels (A-C) show experimentally derived angle-resolved energy loss distributions represented as heat maps for three values of the incidence polar angle θ i indicated with red numbers on the polar axes. The energy loss is the fraction of the incidence kinetic energy of the projectile E i and the kinetic energy after its collision with the surface E s . Panels (D-F) show theoretically predicted distributions derived with the HDNN-PES of this work. In Fig. 6 we also show the distributions obtained when using the REBO-EMFT PES from Ref. 8-see panels (G-I).
The total observed scattering flux decreases rapidly as the incidence angle approaches the normal-this occurs for two reasons. First, the normal component of H/D kinetic energy is more effective in promoting passage over the barrier to chemisorption 8 . Thus, smaller incidence angles produce more sticking. Secondly, the experiment can only observe scattered atoms within a plane defined by the direction of the atomic beam and the normal to the surface. Changing the incidence angle affects the fraction of atoms scattered within that plane. The drop in scattering flux caused by the reduction of the incidence angle is indicated quantitatively by the multiplying factors on the panels. Clearly, the HDNN-PES predictions are in better agreement with experiment than those of the REBO-EMFT PES-see Fig. 6.
Both H and D scattering from graphene exhibit two distinct energy loss channels: a quasi-elastic and a high energy loss channel. The quasi-elastic channel comes from trajectories that fail to cross the barrier to chemisorption, whereas the high energy loss channel arises from trajectories that passed through the chemisorption well forming a transient C-H bond and subsequently returned to the gas phase 8 . The relative intensities of these two channels are also sensitive to incidence angle. The experiment shows that at large incidence angles-see Figs. 6 and 7 panels (A)only quasi-elastic scattering is seen. At small incidence anglespanels (C)-transient chemical bond formation dominates and at intermediate incidence angles-panels (B)-both channels con-  Table 2 in the SI.
tribute to the scattering signal. The angle-resolved energy loss distributions obtained with the HDNN-PES- Fig. 6 panels (D-F)-capture these experimental observations qualitatively better than those obtained with the REBO-EMFT PES- Fig. 6 panels (G-I).
The influence of isotopic substitution on the energy loss spectra can serve as an additional test to validate the accuracy of the HDNN-PES. Comparing the upper panels of Figs. 6 and 7 shows that the experimentally observed branching into the high energyloss channel is somewhat smaller for D than for H under the same incidence conditions. Classical trajectories carried out on the HDNN-PES describe this isotope effect well. Even subtle difference in the angle-resolved energy loss distributions seen in experiment are captured in the trajectory calculations. Compare for example, panels C (experiment) and F (MD with HDNN-PES) of Figs. 6 and 7.
The quality of the results can be more clearly seen in angleintegrated energy loss distributions shown in Fig. 8. Here, the H and D energy loss distributions have in each case been normalized to the integrated scattering intensity of the θ i = 59.5 • distributions. The integrated scattering intensity drops off somewhat  Table 2 in the SI. Fig. 8 Comparing Theory with experiment: Angle integrated energy loss spectra. All incidence conditions are the same as in Fig. 6 too rapidly with decreasing incidence angle, again likely reflecting the influence of out-of-plane scattering. The theoretically predicted energy loss in the quasi-elastic channel is somewhat larger than seen in experiment and the error is larger for H than for D. This might be a quantum effect allowing the H atom to sample the PES closer to the chemisorption barrier producing more inelasticity.
The theoretically predicted D-atom energy loss for the transient bond-forming channel matches experiment remarkably well-see in particular panel D-but the H energy loss is smaller than experimental observation. We note that this is consistent with possible electronically non-adiabatic dynamics where the H atom energy loss is larger than that of the D-atom due to its higher velocity.

Conclusions
We have developed a high-dimensional neural network potential energy surface for H-and D-atoms interacting with a free standing graphene sheet. The potential reproduces a large set of DFT-GGA electronic structure data with high accuracy and is sufficiently efficient to be used in large-scale molecular dynamics simulations. By computing several hundred thousand classical trajectories we demonstrated the utility of the PES by simulating angle-and energy-resolved H-and D-atom scattering experiments similar to those recently published 8 .
The theoretical distributions are remarkably close to those seen in experiment. They accurately capture the branching between a quasi-elastic channel that samples only the physisorption well and a high-energy-loss channel that results from trajectories that traverse the chemisorption well. The simulations also capture subtle differences between H and D scattering seen in experiment that appear as broadening of the angular distribution at specific values of θ i . These results suggest that for scattering at 1.92 eV, neglecting the Pt substrate in the model of scattering dynamics does not introduce large errors.
We do, however, still see systematic differences between experiment and theory. These may be due to failure of the classical or Born-Oppenheimer approximations or both. It is, of course, possible that improved electronic structure data as well as a proper inclusion of the influence of the Pt substrate could explain the remaining discrepancies between experiment and theory. While the present work is only a first step, it demonstrates a crucial milestone toward developing a first principles quality PES that includes the influence of the metal substrate.

Conflicts of interest
There are no conflicts to declare.  Table 2 Showing the total number of simulated trajectories N total , scattered trajectories within detection limit compared to experimental setup N 3 • , the normal component of incidence energy and sticking probability S 0 for H/D at incidence polar angle θ i .

Trajectories visualized
To get a feeling how the scattering happens it is often useful to visualize trajectory as shown in Fig. 5 on the main body of the paper. Fig. 3 provides a side view additional to the top view shown in Fig. 5 of the main text. We also created several movies, where one can follow the trajectories in more detail. We added a top and a side view for the AIMD and HDNNP trajectory. Furthermore, we added three movies showing an example of the quasi-elastic energy loss channel (fast channel), the high energy loss channel (slow channel) and adsorption of the projectile on the surface, which is not possible to see in the experiment. All the movies included in the SI are created using OVITO version 2.9.0.