An experimentally validated neural-network potential energy surface for H atoms on free-standing graphene in full dimensionality
Sebastian Wille, Hongyan Jiang, Oliver Bünermann, Alec M. Wodtke, Jörg Behler, Alexander Kandratsenka
11 INTRODUCTION
An experimentally validated neural-network potentialenergy surface for H atoms on free-standing graphenein full dimensionality
Sebastian Wille a , c , Hongyan Jiang a , Oliver Bünermann a , b , d , Alec M. Wodtke a , b , d , JörgBehler c , d , and Alexander Kandratsenka ∗ a We present a first principles-quality potential energy surface (PES) describing the inter-atomicforces for hydrogen atoms interacting with free-standing graphene. The PES is a high-dimensionalneural network potential that has been parameterized to 75 945 data points computed withdensity-functional theory employing the PBE-D2 functional. Improving over a previously publishedPES [Jiang et al., Science , 2019, , 379], this neural network exhibits a realistic physisorptionwell and achieves a 10-fold reduction in the RMS fitting error, which is 0.6 meV/atom. We usedthis PES to calculate about 1.5 million classical trajectories with carefully selected initial condi-tions to allow for direct comparison to results of H- and D-atom scattering experiments performedat incidence translational energy of 1.9 eV and a surface temperature of 300 K. The theoreticallypredicted 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 dis-crepancies between experiment and theory are likely due to the influence of the Pt substrate onlypresent in the experiment.
H-atom chemisorption to graphene is relevant to hydrogen stor-age , the catalytic formation of molecular hydrogen in the in-terstellar medium and — because hydrogenation of graphenecan induce a band-gap — two-dimensional semiconductor mate-rials . Recently, a full-dimensional PES was reported using firstprinciples energies obtained from Embedded Mean-field Theory(EMFT) to parameterize a second generation Reactive Empir-ical Bond Order (REBO) function . Using classical and semi-classical dynamics calculations, qualitative agreement was ob-tained with H-atom scattering experiments carried out at inci-dence translational energies E i of 1.9 and 1 eV. Furthermore, thetrajectories provided an atomic scale movie at the femtosecondtime scale showing the formation of a covalent chemical bond .The sticking probability could also be calculated using the REBO-EMFT PES and compared well with experiment at E i = eV. Thissuggests that the REBO-EMFT PES is the best available represen-tation of interatomic forces in the H/graphene system.Despite the progress made in that work, two problems re-mained. First, the EMFT data was derived from a model offree-standing graphene, while the experiment was carried out ongraphene that had been grown on Pt(111) . To account for theinfluence of Pt in the simulations, a Lennard-Jones (LJ) interac-tion model with the Pt substrate was included for each atom inthe graphene layer. This improved agreement with experiment, a Department of Dynamics at Surfaces, Max-Planck-Institute for Biophysical Chemistry,Am Faßberg 11, 37077 Göttingen, Germany. E-mail: [email protected] b Institute for Physical Chemistry, Georg-August University of Göttingen, Tammann-straße 6, 37077 Göttingen, Germany c Institute for Physical Chemistry, Theoretische Chemie, Georg-August University of Göt-tingen, Tammannstraße 6, 37077 Göttingen, Germany. d International Center for Advanced Studies of Energy Conversion, Georg-August Uni-versity of Göttingen, Tammannstraße 6, 37077 Göttingen, Germany. suggesting the influence of the substrate may be important. Un-fortunately, it is unclear how to reparameterize the analyticalREBO-EMFT PES from first-principles energies that include thePt substrate. Therefore, the role of the substrate remains uncer-tain. The second problem concerns the fitting error (7 meV/atom)as the REBO function is not flexible enough to closely reproduceelectronic structure data . This complicates the evaluation of thequality of different electronic structure methods, since the fittingerror can easily be larger than the energy differences betweenthe methods being compared. Clearly, a full-dimensional first-principles PES where fitting errors are small and where the roleof the Pt substrate is included would be a significantly better ap-proach to this problem. For both of these problems a solution isoffered by atomistic potentials employing machine learning (ML)methods.In recent years, ML potentials have become a promising newapproach to construct PESs of first-principles quality . Theyhave a uniquely flexible functional form that allows the accu-rate reproduction of reference data sets obtained in electronicstructure calculations, without sacrificing the efficiency neededwhen they are repetitively evaluated in large-scale molecular dy-namics (MD) simulations. ML potentials have been developedfor many systems. These include free-standing and multi-layergraphene as well as graphite and amorphous carbon ,which are closely related to this work. A frequently used typeof machine learning potential suitable for large condensed sys-tems is the high-dimensional neural network potential (HDNN-PES) method proposed by Behler and Parrinello in 2007 .In this paper, we present the first HDNN-PES for H atomsinteracting with free standing graphene, which we validateagainst data obtained from H and D scattering experiments usinggraphene grown on Pt(111), experiments that are similar to thoserecently reported elsewhere. Compared to the REBO-EMFT a r X i v : . [ phy s i c s . c o m p - ph ] J u l COMPUTATIONAL METHODS
PES , we achieve substantially reduced fitting errors without sac-rificing computational performance. Using molecular dynamics,we show that experimentally obtained H/D-atom energy loss andangular distributions are faithfully reproduced. We demonstratethe improvement represented by the HDNN-PES by comparingthe new results to MD simulations done with the previously re-ported REBO-EMFT PES . The remaining deviations between ex-periment and theory likely reflect the absence of the Pt substratein our simulations; however, the influence of the substrate on thescattering distributions appears to be relatively small. The experimental apparatus has been described in detail in Ref.18. H/D-atoms are generated by photodissociating a supersonicmolecular beam of hydrogen/deuterium iodide with a KrF ex-cimer laser producing atoms with incidence energy of 1.92 eV.A small fraction of these atoms passes through two differentialpumping stages, enter the ultra-high vacuum chamber and col-lide with the graphene sample grown in situ on a Pt(111) sub-strate. The sample is held on a six-axis manipulator, allowingvariation of the incidence angle θ i . Recoiling atoms are excited toa long lived Rydberg state ( n = ) by two laser pulses at 121.5 nmand 365 nm via a two-step excitation. These neutral atoms travel25 cm in a field-free region and pass a detector aperture beforethey are field-ionized and detected by a multi-channel plate de-tector. The arrival time is recorded by a multi-channel scalar.The rotatable detector allows data to be recorded at various scat-tering angles θ s . The graphene sample is epitaxially grown ona clean Pt(111) substrate by dosing ethylene (partial pressure × − mbar) at ◦ C for 15 minutes.
High-dimensional neural network potentials (HDNN-PESs) have been the first type of ML potential enabling the simulationsof large condensed systems. In this approach, the total potentialenergy E tot of the system is constructed as a sum of atomic energycontributions, E tot = N atoms ∑ µ = E µ , (1)depending on the local chemical environment defined by a cutoffradius R c , typically in the range between 6 and 10 Å. The positionsof all neighboring atoms inside the cutoff spheres are describedby sets of atom-centered many-body symmetry functions . Theresulting vector of symmetry function values for each atom rep-resents a structural fingerprint that is used as input for an atomicneural network yielding atomic energy contribution E µ into thetotal energy (1). The functional forms of the symmetry functionsensure the necessary invariance of PES with respect to transla-tions and rotations of the system as well as permutations of likeatoms. The atomic neural networks are feed-forward neural net-works and contain a large number of weight parameters, whichserve as fitting parameters for the HDNN-PES. Each element inthe system is modeled by a separate atomic neural network witha specific architecture and values of the weight parameters cal- culated once for each atom of the respective element in the sys-tem. The values of these parameters are determined in an itera-tive optimization process by minimizing the errors of the energiesand forces for a reference data set of representative structuresobtained from electronic structure calculations, typically density-functional theory. Additional structures that may be required inthe reference set in regions of the PES that are not well sam-pled can be suggested by an automatic procedure employing acommittee of HDNN-PESs and a comparison of predicted ener-gies and forces leading to a self-consistent and unbiased gener-ation of the data set. Once a set of weight parameters accuratelyreproducing the reference data has been found, the PES under-goes a series of careful validation steps . Then, the HDNN-PESis ready for applications. For all the details about the HDNN-PES method, the determination of the weight parameters and itsvalidation procedures, the interested reader is referred to severalrecent reviews The Vienna ab initio simulation package (VASP) version5.3.5 has been employed for the reference electronic struc-ture calculations to generate the training set for the HDNN-PES. Density functional theory (DFT) at the generalized gradi-ent approximation (GGA) level of theory using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional with a plane-wave basis has been used in combination with Grimme D2 vander Waals (vdW) corrections . We made use of the ProjectorAugmented Wave (PAW) approach to model the core and va-lence electron interactions. The kinetic energy cutoff has beenset to 400 eV. The Monkhorst-Pack scheme with a × × Γ -centered k -point mesh for the × surface cell has been usedto sample the surface Brillouin zone. With two atoms per primi-tive unit cell, the slab consists of 24 carbon atoms in total and is8.55 Å × using the default value of0.2 eV as the smearing parameter. The threshold for the changein energy between iteration steps when relaxing the electronicdegrees of freedom has been − eV. The iterative procedure described in detail elsewhere was usedto generate the reference data. Briefly, step by step new DFT en-ergies and forces are added for geometries where the HDNN-PESfit does not show the desired accuracy or covers the full config-urational space. These geometries are identified by comparingthe results of several generated HDNN-PESs with differing net-work structures. The reference data set is then extended with theadditional data until convergence is reached.The initial data set consists of energies and forces obtained RESULTS 3.4 Construction of the Neural Network Potential
Fig. 1 Primitive cell containing two C-atoms used to create the × graphene slab. Important high-symmetry sites are indicated by smallwhite balls. This figure and the ones showing surface structures ofgraphene are created using OVITO version 2.9.0 from DFT calculations for about × reference configurations,which were picked up from: (i) ab initio molecular dynamics(AIMD) trajectories simulating H-atom scattering from grapheneat incidence energy of 1.9 eV and incidence angles of 34 ◦ and52 ◦ at surface temperatures of 300 K and 600 K; (ii) geometriesclose to the minimum energy path to adsorption, where the H-atom was put at the lateral position of the C-atom and the z -coordinates were varied over a range of − . Å ≤ H z ≤ . Å and − . Å ≤ C z ≤ . Å, respectively, with 0.025 Å step and withoutstructures with r CH < . Å, whereas the remaining C-atoms werekept at their equilibrium positions; (iii) graphene geometries cho-sen randomly from an AIMD trajectory thermalized at 300 K withH-atom over the surface. The position of the H-atom is chosenrandomly over the whole simulation cell where the z -coordinateranges from 1 to 6 Å. The configuration with a C–H distance of 6 Åand a fully relaxed graphene surface was used as the asymptoticenergy reference. This structure is our energy zero point.The HDNN-PESs fitted to the initial reference data set werethen improved on the set of about . × configurationsobtained from MD simulations of H-atom scattering from agraphene sheet at incidence energy of 1.9 eV in the wide rangeof incidence angles (from 0 ◦ to 90 ◦ in 10 ◦ step) as well as at in-cidence energy of 6 eV and normal incidence angle with surfacetemperatures of 0 K and 600 K starting over high-symmetry sitesshown in Fig. 1. Moreover, we also trained the HDNN-PESs onthe configurations taken from equilibrium MD simulations of thegraphene surface in the wide range of temperatures from 0 K to2000 K. The high-temperature configurations are useful, since thesurface can be heated locally in the neighborhood of the collisionsite.In total, the final HDNN-PES was trained on the reference dataset of 75 945 configurations. The HDNN-PES has been constructed using the RuNNer code. The atomic neural network’s architecture consists of twohidden layers with 15 neurons per layer providing the energies both for hydrogen and carbon atoms. The parameters of symme-try functions are listed in Table 1 of SI. The symmetry functionvalues have been rescaled to the range from 0 to 1. Randomlyselected 90% of the reference data were used to train the NN,whereas the remaining 10% were used as an independent test setto validate the fit and to check for overfitting. The NN weightparameters were determined from the DFT energies and forcesemploying the adaptive global extended Kalman filter . The ini-tial values of the weight parameters have been chosen randomlyin the interval from − to 1. For the weights, a precondition-ing scheme was applied to reduce the initial root-mean-squareerror (RMSE) . The training data in each of the 200 iterations(epochs) of the fit were presented in a random order to reducethe probability of getting trapped in local minima. MD simulations of H-atom scattering from graphene were per-formed using MDT2 code developed to study atomic scatter-ing from various surfaces . The RuNNer subroutines imple-menting the HDNN-PES providing the energies and forces wereintegrated into the MDT2 code. All the results shown in this pa-per have been obtained from MD simulations carried out usingthe RuNNer-MDT2 interface.MD simulations of the H/D scattering from graphene have beencarried out in the NVE ensemble using the standard velocity Ver-let algorithm with a time step of 0.1 fs. The trajectorieswere started with an H-atom randomly put at height of 3.5 Åover the surface and were terminated either when the scatteredatom distance from the surface became larger than 3.6 Å or whenthe trajectory duration exceeded 200 fs. The initial geometry forthe graphene layer was randomly selected from 1000 configura-tions obtained after the equilibration of the surface at 300 K withAndersen thermostat . Those configurations were extractedfrom a single 100 ps-long trajectory with a period of 100 fs. Inthe experiment, graphene is not a single crystal but a composi-tion of two equally abundant orientational domains. The twodomains have a rotational distribution with a Gaussian width of5 ◦ and they are rotated by 27 ◦ with respect to each other . Thisresults in a H-atom velocity vector that is oriented symmetricallywith respect to the two domains. To achieve scattering condi-tions comparable to experiment, the simulations have been car-ried out by averaging over two domains with incidence azimuth φ i = ± . ◦ , where zero for azimuth angle is aligned with a C=Cbond in graphene. In total, ∼ Fig. 2 shows a two dimensional cut through the convergedHDNN-PES developed in this work reflecting structures near theminimum energy path to chemisorption, where the H-atom ap-proaches directly above a C-atom. A physisorption well can beseen at large H z and a deeper chemisorption well at small H z with C z ≈ . Å. The minimum energy path to chemisorption in- RESULTS
Fig. 2 A cut through the HDNN-PES in the vicinity of the minimumenergy path to chemisorption. The H atom is constrained to lie directlyabove a C-atom. H z and C z indicate the distance of H and C, respec-tively, from the plane of the graphene sheet. The physisorption (+) andchemisorption (+) wells have depths of 9 and 657 meV, respectively. Thebarrier to chemisorption ( × ) has a height of 172 meV. volves both degrees of freedom, demonstrating that the C-atom ispartially re-hybridized from sp to sp at the transition state.The depth of the physisorption well of the HDNN-PES has adepth of 9 meV. This compares well with the DFT energy calcu-lated at this configuration 22 meV. The global physisorption min-imum has the H-atom centered over the 6-ring. Here, the HDNN-PES gives the well depth of 11 meV at an H z = . Å. This com-pares reasonably well with the experimentally determined ph-ysisorption well depth (40 meV) and a correlated, counterpoisecorrected wave function calculation of the hydrogen-coronenesystem, which found the minimum at H z = . Å . The previousREBO-EMFT PES had no physisorption well.The chemisorption well depth of the HDNN-PES (657 meV) alsocompares well with DFT (676 meV) but is deeper than that ofthe REBO-EMFT PES (610 meV). Furthermore, the DFT barrier(160 meV) is reproduced well by the HDNN-PES (172 meV) but islower than that of the REBO-EMFT PES (260 meV).The improved quality of the HDNN-PES in comparison to theREBO-EMFT PES is due both to the use of a dispersion correctedfunctional as well as to reduced fitting error. The RMSE fittingerror of the REBO function to the DFT training data was re-ported to be ≈ meV/atom ; furthermore, the REBO functioncannot represent a physisorption well. The flexibility of the neu-ral network—the RMSE for the HDNN-PES is ≈ . meV/atom forenergies in training and test set and ≈ meV/Å for forces intraining and test set, respectively—easily leads to a physically re-alistic physisorption well and a more accurate representation ofthe DFT energies and forces. Fig. 3 shows the fitting error to theDFT energies graphically. While the errors are not randomly dis-tributed, there is no reason to suspect systematic problems withthe PES over the energy range of 10 eV.Figs. 4 and 5 show perhaps in the most impressive way thequality of the NN fitting. Here, two classical trajectories arerepresented, one performed with the HDNN-PES and one with Fig. 3 Fitting error of E HDNN-PES to E DFT . The upper panel shows thecomparison of the two energies and lower panel shows the signed error.DFT energy scale has its zero at configuration corresponding to a relaxedgraphene sheet at T = K with an H atom 6 Å away from the plane ofthe graphene. Δ E p o t / e V ΔFTHΔNNP
Fig. 4 Potential energies from AIMD and HDNN-PES trajectories with E i = . eV, θ i = ◦ and φ i = ◦ . The two trajectories were launchedwith identical initial conditions and both traverse the chemisorption wellbefore returning to the gas phase after a single bounce. The distance ofclosest approach is below r CH = . Å. Movies of the two trajectories canbe found in the SI.
AIMD. The trajectories correspond to the same initial conditionsand are typical of those that will be compared to experiment be-low. Fig. 4 shows the potential energy change along the trajec-tory while Fig. 5 shows the H atom motion in the trajectory ina perspective drawing. The trajectories obtained with these twoapproaches are nearly identical—note that if the fitting were per-fect they would be identical. We now turn to the question: Howwell is the experiment reproduced by classical MD on the newHDNN-PES.Figs. 6 and 7 show comparisons between experiment and the-ory for H and D scattering from graphene, respectively. In both RESULTS
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 coloredball. The divergence between the two trajectories is due to residual errorin the NN fit to the DFT data. A "side view" of the trajectories can befound in the SI. figures, panels (A–C) show experimentally derived angle-resolvedenergy loss distributions represented as heat maps for three val-ues of the incidence polar angle θ i indicated with red numbers onthe polar axes. The energy loss is the fraction of the incidencekinetic energy of the projectile E i and the kinetic energy after itscollision with the surface E s . Panels (D–F) show theoretically pre-dicted distributions derived with the HDNN-PES of this work. InFig. 6 we also show the distributions obtained when using theREBO-EMFT PES from Ref. 8—see panels (G–I).The total observed scattering flux decreases rapidly as the inci-dence angle approaches the normal—this occurs for two reasons.First, the normal component of H/D kinetic energy is more ef-fective in promoting passage over the barrier to chemisorption .Thus, smaller incidence angles produce more sticking. Secondly,the experiment can only observe scattered atoms within a planedefined by the direction of the atomic beam and the normal tothe surface. Changing the incidence angle affects the fraction ofatoms scattered within that plane. The drop in scattering fluxcaused by the reduction of the incidence angle is indicated quan-titatively by the multiplying factors on the panels. Clearly, theHDNN-PES predictions are in better agreement with experimentthan those of the REBO-EMFT PES—see Fig. 6.Both H and D scattering from graphene exhibit two distinct en-ergy loss channels: a quasi-elastic and a high energy loss channel.The quasi-elastic channel comes from trajectories that fail to crossthe barrier to chemisorption, whereas the high energy loss chan-nel arises from trajectories that passed through the chemisorptionwell forming a transient C–H bond and subsequently returnedto the gas phase . The relative intensities of these two chan-nels are also sensitive to incidence angle. The experiment showsthat at large incidence angles—see Figs. 6 and 7 panels (A)—only quasi-elastic scattering is seen. At small incidence angles—panels (C)—transient chemical bond formation dominates and atintermediate incidence angles—panels (B)—both channels con- Fig. 6 Comparing Theory with Experiment for H-scattering fromgraphene at incidence kinetic energy E i = . eV. The energy loss is thefraction of E i and the kinetic energy of the H-atom after its collisionwith the surface E s . Experimental distribution are shown in panels (A–C) along with theoretical distribution found from MD simulations usingthe HDNN-PES (D–F) and the REBO-EMFT PES from Ref. 8 (G–H).The red labeled ticks indicate both the incidence and specular scatteringangles. The integrated signals of panel A, D and G are normalized to 1.The number of trajectories used for the plots are shown in Table 2 in theSI. tribute to the scattering signal. The angle-resolved energy lossdistributions obtained with the HDNN-PES—Fig. 6 panels (D–F)—capture these experimental observations qualitatively betterthan those obtained with the REBO-EMFT PES—Fig. 6 panels (G–I).The influence of isotopic substitution on the energy loss spec-tra can serve as an additional test to validate the accuracy of theHDNN-PES. Comparing the upper panels of Figs. 6 and 7 showsthat the experimentally observed branching into the high energy–loss channel is somewhat smaller for D than for H under thesame incidence conditions. Classical trajectories carried out onthe HDNN-PES describe this isotope effect well. Even subtle dif-ference in the angle-resolved energy loss distributions seen in ex-periment are captured in the trajectory calculations. Compare forexample, panels C (experiment) and F (MD with HDNN-PES) ofFigs. 6 and 7.The quality of the results can be more clearly seen in angle-integrated energy loss distributions shown in Fig. 8. Here, the Hand D energy loss distributions have in each case been normal-ized to the integrated scattering intensity of the θ i = . ◦ distri-butions. The integrated scattering intensity drops off somewhat OTES AND REFERENCES
Fig. 7 Comparing Theory with Experiment for D-scattering fromgraphene at E i = . eV. Experimental distribution are shown in panels(A–C) along with theoretical distribution found from MD simulations us-ing the HDNN-PES (D–F). The number of trajectories used for the plotsare shown in Table 2 in the SI.Fig. 8 Comparing Theory with experiment: Angle integrated energy lossspectra. All incidence conditions are the same as in Fig. 6 too rapidly with decreasing incidence angle, again likely reflect-ing the influence of out-of-plane scattering. The theoretically pre-dicted energy loss in the quasi-elastic channel is somewhat largerthan seen in experiment and the error is larger for H than forD. This might be a quantum effect allowing the H atom to sam-ple the PES closer to the chemisorption barrier producing moreinelasticity. The theoretically predicted D-atom energy loss for the transientbond-forming channel matches experiment remarkably well—seein particular panel D—but the H energy loss is smaller than exper-imental observation. We note that this is consistent with possibleelectronically non-adiabatic dynamics where the H atom energyloss is larger than that of the D-atom due to its higher velocity. We have developed a high-dimensional neural network poten-tial energy surface for H- and D-atoms interacting with a freestanding graphene sheet. The potential reproduces a large setof DFT-GGA electronic structure data with high accuracy and issufficiently efficient to be used in large-scale molecular dynam-ics simulations. By computing several hundred thousand classicaltrajectories we demonstrated the utility of the PES by simulatingangle- and energy-resolved H- and D-atom scattering experimentssimilar to those recently published .The theoretical distributions are remarkably close to those seenin experiment. They accurately capture the branching betweena quasi-elastic channel that samples only the physisorption welland a high-energy-loss channel that results from trajectories thattraverse the chemisorption well. The simulations also capturesubtle differences between H and D scattering seen in experimentthat appear as broadening of the angular distribution at specificvalues of θ i . These results suggest that for scattering at 1.92 eV,neglecting the Pt substrate in the model of scattering dynamicsdoes not introduce large errors.We do, however, still see systematic differences between exper-iment and theory. These may be due to failure of the classical orBorn–Oppenheimer approximations or both. It is, of course, pos-sible that improved electronic structure data as well as a properinclusion of the influence of the Pt substrate could explain theremaining discrepancies between experiment and theory. Whilethe present work is only a first step, it demonstrates a crucialmilestone toward developing a first principles quality PES thatincludes the influence of the metal substrate. Conflicts of interest
There are no conflicts to declare.
Acknowledgements
We thank the DFG for financial support via the SFB 1073 projectsA04 and C03 (Project No. 217133147). We gratefully ac-knowledge the funding of this project by computing time pro-vided by the DFG project INST186/1294-1 FUGG (Project No.405832858).
Notes and references
NATURE , 2001, , 353–358.2 L. Hornekaer, A. Baurichter, V. Petrunin, D. Field and A. Luntz,
SCIENCE , 2003, , 1943–1946.3 R. Balog, B. Jorgensen, L. Nilsson, M. Andersen, E. Rienks,M. Bianchi, M. Fanetti, E. Laegsgaard, A. Baraldi, S. Lizzit,Z. Sljivancanin, F. Besenbacher, B. Hammer, T. G. Pedersen,P. Hofmann and L. Hornekaer,
NATURE MATERIALS , 2010, ,315–319. OTES AND REFERENCES NOTES AND REFERENCES
JOURNAL OF CHEMICAL THEORY AND COMPUTA-TION , 2015, , 568–580.5 F. Ding, T. Tsuchiya, F. R. Manby and T. F. Miller, III, JOURNALOF CHEMICAL THEORY AND COMPUTATION , 2017, , 4216–4227.6 F. Ding, F. R. Manby and T. F. Miller, III, JOURNAL OF CHEM-ICAL THEORY AND COMPUTATION , 2017, , 1605–1615.7 D. Brenner, O. Shenderova, J. Harrison, S. Stuart, B. Niand S. Sinnott, JOURNAL OF PHYSICS-CONDENSED MATTER ,2002, , 783–802.8 H. Jiang, M. Kammler, F. Ding, Y. Dorenkamp, F. R. Manby,A. M. Wodtke, T. F. Miller, III, A. Kandratsenka and O. Büner-mann, SCIENCE , 2019, , 379–382.9 J. Behler,
J. Chem. Phys. , 2016, , 170901.10 F. NoÃl’, A. Tkatchenko, K.-R. Müller and C. Clementi,
Ann.Rev. Phys. Chem. , 2020, , 361–390.11 P. Rowe, G. Csányi, D. Alfè and A. Michaelides, Phys. Rev. B ,2018, , 054303.12 M. Wen and E. B. Tadmor, Phys. Rev. B , 2019, , 195419.13 R. Z. Khaliullin, H. Eshet, T. D. Kühne, J. Behler and M. Par-rinello,
Phys. Rev. B , 2010, , 100103.14 R. Z. Khaliullin, H. Eshet, T. D. Kühne, J. Behler and M. Par-rinello, Nature Materials , 2011, , 693–697.15 V. L. Deringer and G. Csanyi, Phys. Rev. B , 2017, , 094203.16 J. Behler and M. Parrinello, Phys. Rev. Lett. , 2007, , 146401.17 H. Jiang, X. Tao, M. Kammler, F. Ding, A. M. Wodtke,A. Kandratsenka, T. F. Miller III and O. Bünermann, arXiv:2007.03372 [physics.chem-ph] , 2020.18 O. Bünermann, H. Jiang, Y. Dorenkamp, D. J. Auerbach andA. M. Wodtke, Rev. Sci. Inst. , 2018, , 094101.19 J. Behler, J. Chem. Phys. , 2011, , 074106.20 N. Artrith and J. Behler,
Phys. Rev. B , 2012, , 045439.21 J. Behler, Int. J. Quantum Chem. , 2015, , 1032–1050.22 J. Behler,
Angew. Chem. Int. Ed. , 2017, , 12828.23 J. Behler, J. Phys.: Condens. Matter , 2014, , 183001.24 G. Kresse and J. Furthmüller, Phys. Rev. B , 1996, , 11169.25 G. Kresse, J. Non-Cryst. Solids , 1995, ,222â˘A¸S229.26 G. Kresse and J. Furthmüller,
Comput. Mater. Sci. , 1996, , 15â˘A¸S50.27 G. Kresse and D. Joubert, Phys. Rev. B , 1999, , 1758–1775.28 J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. , 1996, , 3865–3868.29 S. Grimme, Journal of computational chemistry , 2006,
27 15 ,1787–99.30 P. E. Blöchl,
Phys. Rev. B , 1994, , 17953–17979.31 H. J. Monkhorst and J. D. Pack, Phys. Rev. B , 1976, , 5188–5192.32 P. E. Blöchl, O. Jepsen and O. K. Andersen, Phys. Rev. B , 1994, , 16223–16233.33 A. Stukowski, MODELLING AND SIMULATION IN MATERIALSSCIENCE AND ENGINEERING , 2010, , 015012.34 N. Artrith and J. Behler, Phys. Rev. B , 2012, , 045439.35 J. Behler, RuNNer - A Neural Network Code for High-Dimensional Potential-Energy Surfaces , Universität Göttingen2020, .36 T. B. Blank and S. D. Brown,
J. Chemometrics , 1994, , 391–407.37 D. J. Auerbach, S. M. Janke, M. Kammler, S. Kan-dratsenka and S. Wille, Molecular Dynamics Tian Xia 2(MDT2); program for simulating the scattering of atoms andmolecules from a surface (GitHub repository). Available athttps://github.com/akandra/md_tian2 , 2020.38 S. M. Janke, D. J. Auerbach, A. M. Wodtke and A. Kandrat-senka,
JOURNAL OF CHEMICAL PHYSICS , 2015, , 124708.39 O. Bünermann, H. Jiang, Y. Dorenkamp, A. Kandratsenka,S. M. Janke, D. J. Auerbach and A. M. Wodtke,
SCIENCE ,2015, , 1346–1349.40 M. P. Allen and D. J. Tildesley,
Computer Simulation of Liquids ,Clarendon Press, USA, 1989.41 W. C. Swope, H. C. Andersen, P. H. Berens and K. R. Wilson,
J. Chem. Phys. , 1982, , 637.42 D. Frenkel and B. Smit, Understanding Molecular Simulation:From Algorithms to Applications , Academic Press, San Diego,2nd edn, 2002, vol. 1.43 H. C. Andersen,
J. Chem. Phys. , 1980, , 2384.44 E. Ghio, L. Mattera, C. Salvo, F. Tommasini and U. Valbusa, The Journal of chemical physics , 1980, , 556–561.45 M. Bonfanti, R. Martinazzo, G. F. Tantardini and A. Ponti, TheJournal of Physical Chemistry C , 2007, , 5825–5829. MD TRAJECTORIES
An experimentally validated neural-network potentialenergy surface for H atoms on free-standing graphenein full dimensionality † Sebastian Wille a , c , Hongyan Jiang a , Oliver Bünermann a , b , d , Alec M. Wodtke a , b , d , JörgBehler c , d , and Alexander Kandratsenka ∗ a SUPPLEMENTARY INFORMATION1 HDNNP details
The NNs for hydrogen and carbon atoms were constructed usingsymmetry functions of radial type ρ m = ∑ n f c ( r mn ) exp h − η r mn i , and angular type φ m = ∑ k , n = m ( + λ cos θ kmn ) ζ ζ − f c ( r km ) f c ( r mn ) f c ( r kn ) e − η ( r km + r mn + r kn ) centered on atom m . Here r mn denotes the distance betweenatoms m and n , θ kmn is the angle between vectors r mk and r mn ; f c is a cutoff function, which gives zero if its argument is largerthan a and 1 otherwise. Indices in sums run over all the neigh-boring atoms of central atom m . η , λ and ζ are the parametersdefining a symmetry function, their values are listed in Table 1.We use 30 input neurons for an H-atom and 60 per C-atom.The quality of the fit to the DFT data is shown for the energiesin Fig. 3 of the main body of the paper. In this supplementarysection, Fig. 1 shows in turn the correlation of the amplitudes offorces | F DFT | extracted from the DFT data and forces | F NN | ob-tained with HDNNP. Fig. 2 presents the same information in theform of the histogram providing a clearer representation of the fiterrors. The RMSE associated with the forces is small, ∼
90 meV/Å.
The most important performance feature of the HDNN-PES is itsability to accurately calculate the energy and forces for a sys-tem with many degrees of freedom with low computational costs.This is particularly important in the simulations of angle resolvedatomic scattering experiments from surfaces. The experimentaldesign detects only a small fraction of the scattered atoms due toits high angular resolution. Consequently, one needs hundreds ofthousand trajectories to reduce statistical noise to the level of the a Department of Dynamics at Surfaces, Max-Planck-Institute for Biophysical Chemistry,Am Faßberg 11, 37077 Göttingen, Germany. E-mail: [email protected] b Institute for Physical Chemistry, Georg-August University of Göttingen, Tammann-straße 6, 37077 Göttingen, Germany c Institute for Physical Chemistry, Theoretische Chemie, Georg-August University of Göt-tingen, Tammannstraße 6, 37077 Göttingen, Germany. d International Center for Advanced Studies of Energy Conversion, Georg-August Uni-versity of Göttingen, Tammannstraße 6, 37077 Göttingen, Germany. † Electronic Supplementary Information (ESI) available: [details of anysupplementary information available should be included here]. See DOI:10.1039/cXCP00000x/
Table 1 Parameters for the atom-centered symmetry functions for H- andC-atoms. Equations for the types used can be found in section 1. no. η ( a − ) λ ζ Radial1 0.0002 0.0053 0.0134 0.0275 0.0606 0.156Angular7 0.000 1 18 0.000 1 29 0.000 1 410 0.000 1 1611 0.000 -1 112 0.000 -1 213 0.000 -1 414 0.000 -1 1615 0.013 1 116 0.013 1 217 0.013 1 418 0.013 1 1619 0.013 -1 120 0.013 -1 221 0.013 -1 422 0.013 -1 1615 0.156 1 116 0.156 1 217 0.156 1 418 0.156 1 1619 0.156 -1 120 0.156 -1 221 0.156 -1 422 0.156 -1 16experimental data. Table 2 shows the total number of MD trajec-tories calculated at the conditions used in the paper as well as thenumber of trajectories scattered into the cone with the apex angleof ◦ corresponding to the geometry of the experimental setup.When simulating hydrogen scattering from graphene, the tra-jectories were interrupted after a 200 fs simulation time. H atomswhich had not left the surface after this time were considered tohave adsorbed. The sticking probabilities for the incidence condi-tions used in the main paper are collected in Table 2. a r X i v : . [ phy s i c s . c o m p - ph ] J u l TRAJECTORIES VISUALIZED
Fig. 1 Correlation plot of HDNNP | F HDNN-PES | and DFT forces | F DFT | for the training set in blue and test set in red, respectively. The HDNNPrepresents the DFT forces over the whole range quite well or at leastreasonable. −6 −4 −2 0 2ΔFΔ(eVΔ/ΔÅ)0100000200000300000400000500000600000700000800000 C o un t s TrainingΔSetΔmeanTestΔSetΔmeanTrainingΔSetTestΔSet
Fig. 2 Histogram of differences of the forces F DFT − F HDNN-PES for thetraining set in blue and test set in red, respectively. The mean values forthe training and test set are also shown in the corresponding colour.Table 2 Showing the total number of simulated trajectories N total , scat-tered trajectories within detection limit compared to experimental setup N ◦ , the normal component of incidence energy and sticking probability S for H/D at incidence polar angle θ i . θ i N total N ◦ E i cos θ i /eV S H 40 354,911 22,630 1.13 0.2150 291,238 22,533 0.79 0.3959.5 322,691 122,205 0.49 0.22D 43 189,577 7,021 1.03 0.5251 183,837 10,452 0.76 0.6859.5 221,861 73,607 0.49 0.39
Fig. 3 Side view of Fig. 5 of the main text.