Quasiclassical Trajectory Dynamics Study of Atomic Oxygen Collisions on an O-Preadsorbed Graphite (0001) Surface with a New Analytical Potential Energy Surface

A new flexible periodic LEPS potential energy surface (FPLEPS) based on density functional theory data is constructed for the interaction of atomic oxygen with an O-preadsorbed graphite (0001) surface over a C–C bridge. New ingredients were added to the usual expression of the FPLEPS in order to take into account the entrance barriers, molecular orientation, and morphology of the surface. A total of 563 DFT points were used to fit the Eley–Rideal (ER) reaction channel, achieving a root-mean-square deviation of 0.120 eV for energies lower than 1 eV over reactants. A quasiclassical trajectory (QCT) dynamics study has been performed at several initial conditions: collision energies (0.01 ≤ Ecol ≤ 2.0 eV), incident angles (θv = 0°, 45°), and surface temperatures (100 ≤ Tsurf ≤ 900 K). Also quasithermal and hyperthermal (⟨Ecol⟩ = 5.2 eV) conditions were considered. Eley–Rideal reaction and O reflection were the main processes, keeping the formed O2 molecules translationally and internally excited via the ER pr...


I. Introduction
Carbon fiber technology has already been used to replace many spacecraft components. In particular, fiber reinforced materials such as carbon and glass composites have the highest strength and stiffness to weight ratios among engineering materials. Carbon fiber composites, for example, are five times stiffer than steel for the same weight allowing for much lighter structures for the same level of performance. Since the strength of graphite improves with temperature up to about 2,760 K, it is possible that vehicles, which must enter Earth's or Jupiter's atmosphere, or orbit very close to the Sun, may have some structural parts made out of carbon materials as graphite/graphene. Satellites or spacecrafts in low Earth orbits (LEO), at approximately 160-2,000 km above Earth's surface, must travel very fast so that gravity will not pull them back into the atmosphere (i.e., around 7.5 km/s). At these conditions some of the main present gas species (i.e., O/O2) collide with the spacecraft surface [1] at high collision energies (i.e., about 4.5 ± 1.0 eV) producing all kind of heterogeneous reactions (i.e., etching, adsorption, recombination,...). Thus, it is very important for space industry to describe properly the interaction between O/O2 with graphite and other related materials at these extreme conditions. Some experiments [2,3,4] are performed with beams of O( 3 P) atoms with some extent of O2(X 3 Σg -) that impinge over highly ordered/oriented pyrolitic graphite (HOPG) and analyze the chemical reactivity and morphological evolution of the surface. The main conclusion is that the graphite surface becomes functionalized with epoxide groups that can migrate across the surface. Moreover, incoming O atoms react at the surface to produce O2 molecules via an Eley-Rideal (ER) mechanism whereas incoming O2 molecules are only inelastically scattered.
Theoretical works regarding the interaction of O/O2 with some graphitic surfaces (e.g., graphite [5,6,7,8,9,10], graphene [11,12,1314] and carbon nanotubes [7,15,16]) confirm the experimental observations. They point out that impinging O atom becomes chemically adsorbed on a bridge position between two adjacent C atoms of graphite or graphene, forming an epoxide group with an adsorption energy in the range of 0.95 -3.4 eV [5,6,7,8,9,14] for a basal graphite surface depending on the oxygen coverage and the level of the calculations. When the interaction is over a defective edge site on graphite surfaces, the adsorption energy increases considerably (4.5 -7.5 eV) [9,14]. Another important feature of the atomic oxygen adsorption over graphite is the easy diffusion observed (energy barriers of 0.38 eV [5] and 0.36 eV [6] for a coverage of 12.5 % and 3.12 %, respectively), which are lower than the ones obtained for a graphene sheet (0.58 eV for a coverage of 16.7 % [11]), that would imply a noticeable mobility of isolated adsorbed atom. In addition, we recently proposed a microkinetic model for O/O2 over graphite [5]. This model used rate constants obtained from DFT and standard transition state theory for several processes occurring when O/O2 mixtures collide over a clean graphite surface, showing a very low steady-state atomic coverage (q < 0.5%) and also very low atomic recombination coefficients (g) at quasithermal conditions (qv, T = 300-1,000 K).
Hyperthermal collisions at normal incidence of O over graphene [2,13] and graphite [13] have been simulated by using DFT (PBE/DZP) data and the ReaxFF reactive force field, respectively, but without including the thermal bath. The surfaces were functionalized with preadsorbed O atoms or single vacancy defects. These dynamical simulations allow very long integration times (e.g., 150 ps [13]), and also increase the system size (every picosecond a new incoming oxygen atom impinges the surface, increasing thus the system size by one O atom per collision). In these conditions, it was observed that O adsorption (i.e., epoxide and carbonyl formation) and O2 formation via an ER mechanism were the main processes (t < 45 ps) while inelastic collisions were very low. At longer times CO/CO2 formation could also be observed and C sheet would be both broken and deformed. Nevertheless, neither 5 eV oxygen atoms nor 10 eV oxygen molecules are expected to cause sputtering from a pristine unfunctionalized surface [17] at least if no defect sites are previously formed. The absence of CO formation arising from pristine graphene at 3,000 K is also confirmed by DFT tight binding with dispersion (DFTB-D)-based molecular dynamics simulations [14].
The use of an analytical expression for the potential energy surface (PES) allows obtaining easily and efficiently the potential energy and its first order derivatives for performing fast dynamics calculations. In the last decades, the PLEPS (Periodic London-Eyring-Polanyi-Sato) function has been extensively used for studying the dissociative diatomic molecular adsorption mainly over several metallic systems [18,19,20,21,22,23,24,25,26,27,28,29,30] and for the Eley-Rideal (ER) recombination of some particular systems (H/H-Cu(111) [31,32,33,34], H/D-Si(001) [35], H/H-graphite [36,37]). The popularity of PLEPS potential is due to the simplicity of its analytical form. However, the PLEPS model is known to be not flexible enough to describe the intricate structure of high-dimensional PES [38,39,40]. The use of only two adjustable parameters (Sato parameters) to fit a multi-dimensional PES (6D-PES for diatomic/surface systems) explains this major shortcoming. This is particularly true for heavy diatomic molecules interacting with surfaces that exhibit in many cases strongly corrugated PES [41,42]. Thus, although the PLEPS model was one of the reference models in the past, nowadays this method has become unreliable to deal with gas-surface reactions. An extension of the PLEPS model called FPLEPS (Flexible PLEPS) has been recently applied to N2/W(100,110) [40,43,44] and O2/Cu(100) [45], reaching a good level of accuracy. The FPLEPS is flexible enough to be adapted to various atom/diatom surface systems and for the study of various elementary reactions. Moreover, a new FPLEPS was recently provided [46] in parallel to an interpolated Modified Shepard one, for studying the O collisions with a pristine graphite (0001) surface giving dynamical information by means of classical trajectories in good agreement with the available experimental hyperthermal data [2].
The main goal of this work is to provide a new analytical FPLEPS PES for describing the O collisions with an O-preadsorbed graphite (0001) surface by using the previous developed PES for the O-atom interaction with the clean surface [46]. A quasiclassical trajectory (QCT) dynamics study is carried out with this new PES and compared with the available experimental data [2]. This paper is organized as follows: Section 2 provides a brief description about DFT and QCT calculations; Section 3 describes the FPLEPS PES construction and its main characteristics and Section 4 presents the main dynamical results and discussion, including the comparison with experimental data.
Finally, Section 5 gives the summary and conclusions.

A. DFT calculations
Periodic density functional theory (DFT) calculations have been performed using the Vienna ab initio simulation package (VASP) [47,48,49,50]. A detailed description of the DFT calculations was done in previous works [5,46], thus only a few details are given here. The calculations are based on the generalized gradient approximation (GGA) with the revised Perdew-Burke-Ernzerhof functional (RPBE) [51]. The projector-augmented wave (PAW) technique within the frozen core approximation has been used to describe the electron-core interaction [52,53]. An energy cut-off of 550 eV has been used in the plane-wave expansion and the Brillouin zone was integrated by using an (11 ´ 11 ´ 1) kpoint mesh by means of the Monkhorst-Pack method [54]. All the slab calculations used a vacuum around 15 Å, large enough to prevent significant interactions between periodic images. The energy convergence in the electronic self-consistent procedure was maintained below 10 -6 eV for all geometrical calculations performed within the rigid slab model. Due to the important role that spin plays in atomic oxygen description, all calculations were spin-polarized [5,7,9]. However, as the Eley-Rideal reaction does not imply necessarily a change in the spin state because both asymptotes correspond to triplet states (i.e., reactants (O (ad) + O (g) ) and products (O 2(g) + graphite)), this gives further support to the reliability of the present DFT approach. Moreover, our previous DFT study on O/graphite system [5,46] showed a good agreement with available ab initio quantum molecular structure data and previous DFT data. Despite everything, we have carried out additional DFT calculations in order to calibrate the methodology used. Thus, non-periodic DFT calculations have also been performed with either UHF or RHF reference functions, using different basis sets and with a couple of hybrid functionals (PBE and B3PW91), by means of Gaussian 09 [55]. The cluster model used to simulate the graphite surface has been the pyrene molecule (C 16 H 10 , cluster A in Ref. [16]). The results obtained for the adsorption energy, for the ER exothermicity and for the ER barrier are in good agreement with the results reported in this and in our previous work [5], taking into account the lack of periodicity in the cluster model. The differences between our values of adsorption energy and those reported in Refs. [7] and [9] can be rationalized by the fact of using different functionals. However, these DFT uncertainties could be reduced using more accurate ab initio methods (e.g., CASPT2, Coupled Cluster,..), although DFT methods allow an easier treatment to model huge and periodic systems.
According to the previous study [5], molecular channel does not present any particular feature and is essentially repulsive close to the surface. A high endothermicity prevents the molecular dissociation, so the DFT data used to obtain the analytical surface was focused on the ER entrance channel. Since the most stable atomic adsorption site was the bridge one (B) for the atomic PES [46] as well as for DFT data [5], the configurations needed to be included into the PES fitting procedure consider also an oxygen atom (OA) fixed to the DFT calculated bridge minimum for a (3 ´ 3) supercell (i.e., 11.1 % of coverage considering bridge sites). Working with such a big supercell increases considerably the number of atoms to be included in DFT calculations and consequently the computational cost, but it was necessary to avoid interactions between adsorbed atoms in contiguous cells. To ensure a proper description of the interaction of the incoming oxygen atom (OB) with the preadsorbed one, configurations at different f angles (i.e., 0º, 30, 60 and 90º, Fig. 1) were calculated for several impact parameters in the range of 0.0 Å < b < 3.0 Å and considering the incoming atom at ZB values between 1.1 Å and 3.1 Å. Thus, a total of 563 DFT configurations were calculated, whose contour plots are shown in the right panels of Fig. 2.

B. QCT calculations
A dynamical study involving the interaction of atomic oxygen over O-preadsorbed graphite surface was performed by means of the QCT method [56]. Several initial conditions were sampled in order to investigate state-specific, quasithermal and hyperthermal processes. State-specific (E col , Tsurf, θv) calculations were carried out for fixed initial collision energies in the range of 0.01 eV £ E col £ 2.0 eV and at two incident velocity angles (θv = 0° and 45°) for a given surface temperature Tsurf. The atomic incoming velocity angle is defined with respect to the negative Z axis (e.g., 0° for normal incidence) and its projection onto the X-Y plane (i.e., the azimuthal f v angle) was uniformly sampled within the 0° -360° interval.
Initial incoming atomic oxygen position (XB, YB) is randomly selected along the (1 ´ 1) unit cell meanwhile initial ZB position is set to 7.5 Å where the interaction with surface is negligible. The preadsorbed oxygen atom (OA) was fixed initially at the equilibrium geometry of the B1 minimum found for the atomic FPLEPS surface [46] (XA = 0.000 Å, YA = 0.713 Å, ZA = 1.385 Å, Fig. 1). To account for the energy exchanged and dissipated, a GLO model [56,57,58,59] was used, including the surface and ghost particle motion into the Hamilton equations and by using the same optimum parameters, which were justified previously for O/graphite interaction [46]. These parameters were: wi,x = wi,y =10 -3 and wi,z =3.4 ´10 -4 a.u. for i = 1,2,3 (oscillator frequencies), m = 60 amu (effective mass) and gg,x = gg,y = gg,z = 4.0 ´10 -4 a.u. (friction constants). The surface temperatures (Tsurf) studied were in the range of 100 K £ Tsurf £ 1,500 K. The preadsorbed O atom is initially (t = 0) displaced with the same coordinates of the surface particle (Xs, Ys, Zs) to avoid any change in its potential. To make an easier initial selection for adatom coordinates, we put this at its equilibrium geometry (i.e., potential energy V = 0) and then its total energy (i.e., E = 3k B T surf ) was introduced only as kinetic energy (i.e., sampling a Boltzmann distribution with T = 2×Tsurf, [60] and references therein), assuming thus the virial theorem. We checked that this adatom was correctly thermalized together with the GLO bath during the trajectory and before the arrival of the incoming atom. The adsorbed atom remains around the equilibrium position and shows the usual oscillations for the kinetic and potential energies. The total energy is not constant due to the thermal bath dissipation effect.
The qctsurf code developed in our group was used to calculate the trajectories integrating the Hamilton equations of the system (included the GLO) using the Beeman algorithm. The time step used was 5´10 -17 s, which ensures a total energy conservation along the trajectories lower than 1´10 -4 eV in absence of the thermal bath. Total integration time was set to a maximum of 2.5 ps, long enough to obtain one of the possible exit channels defined as followed. Adsorption and reflection of the impinging atom was considered either if the preadsorbed atom remains adsorbed to the surface or if it desorbs, respectively. For atomic adsorption classification, the ZA,B should be lower than 2.2 Å and the number of total rebounds with the surface more than 8, as was also set previously [46] for the atomic case over clean graphite. On the other hand, any of the oxygen atoms were considered reflected if either ZA or ZB were higher than 7.6 Å and the direction of the velocity vector pointed to the vacuum.
The formation of the O2 molecule was also considered through two different channels either if the final molecule remains adsorbed or if it goes to the gas phase (i.e., ER reaction). The oxygen molecule was considered formed when the internuclear distance R was lower than 2.0 Å and dissociated for higher values than 3.0 Å. We also checked the formation of O2 molecules by either a direct (ER) mechanism or by a non-direct (hot atom) one. The criterion for a first tentative classification of some trajectories as ER or hot atom atomic recombination was based in the number of atomic rebounds previous to the O 2 formation. Thus, reactive trajectories with more than two rebounds were classified as non-direct.
Similarly as for atomic adsorption, when Zcm was smaller than 2.2 Å and the number of total rebounds with the surface was more than 8, the event was classified as molecular adsorption.
The minimum number of total trajectories (NT) calculated for each condition was around 10,000, which leads to a standard error of 1 -2% in ER probabilities.

III. Analytical FPLEPS potential energy surface
The complete expression used for describing the 6D potential V(RA,RB) of the two atoms (A, B) interacting with the graphite surface ( Fig. 1) is defined as: (1) where the positions of the two atoms A and B are respectively given by the vectors RA(XA,YA,ZA) and RB(XB,YB,ZB). G1 and G2 are two Gaussian functions introduced to improve the description of the full potential in the entrance valley in the same manner that in a previous work [40]. The terms and represent respectively the Coulomb and exchange integrals for two-body systems and are associated with the atom-surface interaction (i = As, Bs, described previously [46]) and the molecular interaction (i = AB). These integrals have the usual form: with Ri = R for i = AB and Ri = ZA, ZB for i = As, Bs, and where Di, αi and are the well known Morse parameters (i.e., dissociation energy, Morse shape parameter with the introduction of some particular flexibility and the equilibrium distance). The Sato parameters (D As =D Bs , D AB ) are chosen to fit the main features of the ER DFT channel. The periodicity in the surface is achieved by assigning a (X,Y)-dependence to the full set of parameters, concretely by means of Fourier expansions. Moreover, a switching function fi was also added to Ui and Qi only for i = As and Bs (not for i = AB) in order to fit better the atomic adsorption barriers, The dependence in (Xk,Yk) of the full set of parameters in Eq. (4) was ensured by means of a 5 th order Fourier expansion. All the details related with the oxygen-graphite interaction were developed previously [46].
The 2D Gaussian functions (G1 and G2) improve the description of the entrance valley shape in the ER reaction, where Zcm refers to the molecule centre of mass, and b to the impact parameter (Fig. 1a). The other parameters (i.e., Ai, bi, Zcm,i, σi and σcm,i) are f dependent and are related with the height, position and the width of the Gaussian functions, whose optimal values are found by fitting the DFT data for each selected f angle (i.e., 0°, 30°, 60° and 90°), listed in Table I. It was necessary to fit also the Sato parameters for each f angle in order to describe properly the complexity of the ER entrance channel. Fig. 2 shows the good agreement between DFT and FPLEPS contour plots for the four angles used for describing the ER channel. A periodic function (P) has to be introduced to represent the f-dependence of the FPLEPS parameters. This function must accomplish and , which corresponds to the periodicity of the surface when atom A is preadsorbed over a determined bridge site (the B1 site is arbitrarily chosen in this work, cf. Fig. 1b). This function is defined for each parameter Pj (j = 1, 12) as follows, The Bij parameters that solve the equation system are listed in Table II. It is worth noting that the previous strategy works when the preadsorbed atom is located over a particular bridge site (B1) but fails when trying to reproduce the correct behaviour over other bridge sites (B2 or B3). This is due to the different bridge sites are not equivalent with regard to the f value within the framework of two atoms interacting on a graphite surface. Thus, for example the f = 0° orientation over B1 site corresponds to f = 60° and 120° over B2 and B3 sites, respectively (assuming that T1 and T2 sites are fully equivalent). The simultaneous correct description of the three bridge sites requires considering a new Fourier expansion. In order to keep the good symmetry properties, it is necessary to include not only the three bridge sites in the new Fourier expansion but also nine additional sites (i.e., T1, T2, H, T1H1, T1H2, T1H3, T2H1, T2H2 and T2H3 defined in Fig. 1b). The idea is to make also a (X,Y)dependence of the molecular parameters (Sato and Gaussian parameters) via this new Fourier expansion (as for the atomic part but with different symmetry properties related to the two-atom interaction). A Fourier expansion series up to the 11 th order is used, where Cji (i = 0 to 11 and j = 1 to 12) represent the Fourier coefficients, a = b = 2.468 Å is the cell parameter and (u,v) are the crystal coordinates (Fig. 1b), which are related with the orthogonal (X,Y) coordinates by (9) Before using the Fourier expansion for the Sato and Gaussian parameters, it is necessary to know the values of these parameters for the nine other sites that are included in the expansion. It was not necessary to perform additional DFT calculations since over these sites, the atomic adsorption is either weaker than over the bridge sites (T1, T2 sites) or simply not possible (potential is quite repulsive).
Consequently, the value of the different parameters over these sites was not so important for the correct ER channel description. Thus, in order to ensure the surface symmetry requirements, an average value of the parameters for the three bridge sites was used for the nine other sites. Some checks were done in order to know the FPLEPS behaviour under this approximation, but the most important thing was to describe accurately the bridge sites, where the adsorption is much more favourable and where it is very likely that the ER processes happen.
One important aspect was to make sure that the PES was equivalent for any exchange of both atoms. The Finally, as angle f is not defined for b = 0 (i.e., perpendicular configuration), the determination of Pj using Eq. (7) is not possible. This is only true for the Sato parameters since the Gaussian ones are not located in that region of space configuration. Therefore, it was necessary to include a new switching function for the two Sato parameters in the region close to perpendicularity (b ≈ 0). In this region Pj is then defined by the expression, where the top value is used when j = 1 (i.e., for ΔAB) and the lower one in the case of j = 2 (i.e., for ΔAs, Bs,); the Pj term inside the equation is the one obtained in Eq. 10, whereas the Pj on the left part is the final correct value. The switching functions used are, where the numerical values inside the switching functions were chosen in order to give a connection as smooth as possible in the b and R intervals considered. Thus, fb is only different of 1 in the limit b ® 0 Å, otherwise its value is the unity and Pj (j = 1, 2) parameters remain unchanged.
Details about the binding energy of the O atom to the surface over different sites are reported and discussed widely in our previous work [46]. As a summary, only the asymptotic energy values for the resulting FPLEPS PES are listed in Table III.  For the preadsorbed O atom onto B1 site (Fig. 3a), the favoured approaching angles for the incoming atom are 60°, 120°, 240° and 300°. It is important to note that these favoured approaching angles are obtained for ZB = 3.00 Å. Indeed, if the distance to the surface diminishes the approach over the closest neighbour bridges is not preferred at all and the incoming atom tends to get adsorbed over a non-adjacent bridge site (Fig. 3d), as predicted in the previous DFT theoretical work [5]. confirm that this barrier should be lower than 0.04 eV (PBE/cc-pVTZ).

A. State-specific conditions
Probabilities for each of the exit channels were calculated for two incident angles (θv = 0° and 45°) at several initial collision energies between 0.01 eV and 2.0 eV. Results plotted in Fig. 4 show that the reflection of the impinging atom and the ER reaction are the two main processes observed in all the conditions explored.
For normal incidence without thermal bath (Fig. 4a), the ER probability firstly increases with Ecol and latter decreases at high collision energies, without showing threshold energy in agreement with the negligible energy barrier found in the PES. The reflection follows just the contrary behaviour. Atomic adsorption is much less important and becomes open at medium energies while atomic desorption is only produced at very high energies, where the impinging atom destabilizes the preadsorbed one.
When off normal incidence is considered (Fig. 4b), ER reaction probability decreases considerably compared with normal incidence results (Fig. 4a) becoming the atomic reflection more important than the ER reaction. Moreover, desorption processes are now more significant, probably because the parallel component of the kinetic energy of the impinging atom favours the removing of the adatom.
The effects of the surface temperature introduced through the thermal bath (Figs. 4c-f) are small.
Only some differences at high collision energies are observed. Thus, the reflection process is slightly attenuated since the atomic adsorption becomes larger. This fact is due to an appreciable energy transfer from incoming atom to the surface, hence enhancing the atomic trapping.
Since the preadsorbed O atom remains weakly linked to the graphite surface (E ad = 0.68, Table III), the possibility of extracting this adatom from the surface with an incoming atom is easier than when the adatom is strongly adsorbed (E ad = 3.05, 2.35 and 2.8 eV for H over W(100,110) [63], Cu(110) [64] and Pd(111) [65], respectively or E ad The adsorption and reflection probabilities are lower than the ones obtained for a clean graphite (0001) surface [46], because now the preadsorbed oxygen occupies the bridge site (so, there are not any free bridge site inside the (1 ´ 1) unit cell for adsorption) and also due to the existence of other processes (i.e., ER, desorption), that obviously compete with the adsorption and reflection.
The analysis of the initial positions of the impinging atom shows that trajectories leading to ER reaction generally start at positions close to the preadsorbed atom because they allow the interaction between both oxygen atoms to form the O2 molecule, in agreement with PES topology (Fig. 2).
Furthermore, the atomic adsorption process without desorption of the preadsorbed atom is only possible when the R distance is long enough to consider the interaction almost negligible; in this case, the incoming atom is behaving as if it was interacting with a clean surface. The analysis of the final position of these trajectories shows that this atom finally adsorbs over bridge sites located close to the starting position of the trajectories but in the surrounding cells with free bridge sites.
Polar scattering angle (θv') distributions of reflected atoms are plotted in Fig. 5 for three initial collision energies (Ecol = 0.2 eV, 0.5 eV and 1.2 eV), two initial incident angles (θv = 0° and 45°) and different surface temperatures (i.e., without thermal bath, 100 K and 900 K). These distributions can be compared with the collisions over a clean surface [46]. For normal incidence there is a clear shift of the peaks towards higher scattering angles (e.g., ~ 60° for Ecol = 0.20 eV) while for off normal incidence this shift corresponds to smaller angles (e.g., ~ 40° for Ecol = 0.20 eV), indicating clearly that the presence of preadsorbed oxygen atom affects these angular distributions. The analysis of initial position of the impinging atom for normal incidence regarding the final θv' angle shows that atoms reflected at low angles are located initially far from the preadsorbed one while the high angles are mainly originated by atoms reflected close to the adatom. With the clean surface, the entrance barriers were quite similar over the main sites and therefore, the reflection on these barriers did not change a lot the initial incident angle [46]. Now, according to Fig. 2, the position of the barriers depends a lot of the impact parameter. Consequently, these barriers appear to be very corrugated with respect to the incoming atom leading to a broader angular distribution peaked around 60°. Obviously, this new PES topology is clearly created by the preadsorbed atom. Moreover, the higher is the initial Ecol, the lower is this topological effect since the barriers can be overcome more easily. Peaks appear at low angles (10° -20°) for normal incidence and they are a bit shifted at larger ones (15° -40°) for off normal incidence. In general, the effect of surface temperature on the normal incidence distributions (Fig. 6, left panels) is small because the ER process occurs mainly via a direct mechanism and the interaction with the surface is negligible except for the preadsorbed atom that is thermalized with the surface. For off-normal incidence (Fig. 6, right panels) the interaction with the surface is higher and some distributions are double peaked. This fact shows the presence of a non-direct mechanism [67], more evident for lower collision energies at Tsurf = 100 K (Fig. 6d) and practically insignificant when the surface temperature is increased (Fig. 6f). The analysis of trajectories shows that non-direct mechanism produces O2 molecules predominantly formed at low scattering angles (10° -20°), whereas the direct mechanism drives to O2 in a wide range of scattering angles (< 45°).
Moreover, non-direct mechanism is mainly obtained for off-normal incidence whereas for normal one its contribution is practically constant (< 4 %) at all the conditions explored. Concretely, for θv = 45° and Ecol = 0.2 eV, the non-direct formed O2 molecules account for 20% and 13% without thermal bath ( Fig. 6b) and for a Tsurf = 100 K (Fig. 6d), respectively. When Ecol is increased to 1.0 eV at Tsurf = 100 K, the non-direct mechanism becomes a bit less favourable (11%), in agreement with a similar behaviour reported for H recombination over Ni(100) [67]. Obviously, in the present study, the incoming O atom samples initially only a (1 ´ 1) unit cell minimizing the possibility of the non-direct mechanism to occur. Thus, a deeper study is planned for sampling a bigger unit cell at several initial conditions.
The exchange of energy between the impinging atom and the preadsorbed surface (i.e., ∆Ecol = Ecol' -Ecol) was also evaluated at three selected collision energies (Ecol = 0.2, 0.5 and 1.2 eV) with Tsurf = 100 K and 900 K for normal incidence, and shown in Fig. 7. This energy exchange is in general negative (e.g., <ΔEcol> = -0.69 eV for Ecol = 1.2 eV), which means that the final collision energy of the reflected atom is lower than the initial one. The average energy exchange tends to diminish at higher temperatures although with rather broader distributions as could be expected. The behaviour is very similar as the one obtained for collisions over a clean surface [46].
An analysis of the internal states (v',j') and the translational energy (Ecol') of the formed molecules via ER reaction is presented in Fig. 8 at two initial collision energies (Ecol = 0.2 and 1.2 eV) for normal incidence. Only Tsurf = 500 K is considered as surface temperature effect was expected to be almost

B. Quasithermal and hyperthermal conditions
Reaction probabilities at quasithermal conditions (θv, Tsurf = TO = T) for gas temperatures between 300 K and 1,500 K and for θv = 0° and 45°, are presented in Fig. 9. The reflection of the impinging atom is the predominant process being the ER reaction the second one, with a probability of around 10% -30% at the studied temperatures. The rest of the processes are almost negligible, except the atomic reflection with the simultaneous desorption of the preadsorbed atom at high temperatures (i.e., T = 1,500 K). For off normal incidence, the ER probabilities (Fig. 9b) are even lower and hence the reflection probabilities become higher. This behaviour is consistent with the one found for low collision energies (Fig. 4), as Tsurf had a small influence on reaction probabilities. The energy exchange with the surface was small (<ΔEcol> = -2.96×10 -3 , -7.07×10 -3 and -7.94×10 -2 eV for Tsurf = 300, 900 and 1,500 K).
The O2 molecules formed via the ER reaction show angular distributions for normal and off normal atomic incidence with similar trends to those observed for state-specific conditions (Fig. 6). Their final energies (i.e., internal and translational) present similar values to those obtained for Ecol = 0.2 and 1.2 eV (Fig. 8). The mean values of these energies (<Ecol'> = 1.28, 1.52 eV and <Eint ' > = 2.53, 3.28 eV at 300 and 1,500 K, respectively) show that the increase of temperature enhances the total energy excitation of the molecule.
In order to study the possible coverage effect in the final polar angle distribution of reflected O atoms over HOPG in hyperthermal experiments [2] with Tsurf = 503 K and θv =45°, we have calculated a large number of trajectories (1,000,000) simulating exactly the experimental O velocity distribution, whose average translational energy corresponds to 5.2 eV. As HOPG surface consists of many orientations within the topmost layer, we have uniformly sampled the azimuthal velocity angle (0° ≤ fv ≤ 360°).
The calculated distribution of reflected O atoms was split between in-plane and out-of-plane contributions (Fig. 10a), using an angular tolerance of ± 0.5° for the in-plane distribution. The major part of the scattered O atoms were produced out-of-plane (99.6%), hence the in-plane contribution is negligible in front of out-of-plane one to account for the total angular distribution. However, both are normalized to the unity in their maxima to enhance their comparison with experiments, which measured scattered atoms in the same plane as the initial plane formed by the normal to the surface and the O beam line. The QCT in-plane distribution shows a peak at low angles (~35°) far from experiment (~62°). This disagreement cannot be explained by the differences between in-plane and out-of-plane distributions. Nevertheless, the comparison of the QCT distribution for clean graphite matches nicely the experimental one [46], although the experimental HOPG surface coverage [2] is unknown. This seeming discrepancy could be explained by the fact that at hyperthermal conditions the atomic oxygen adsorption over a clean graphite surface is very low (0.15%), as showed by our QCT prediction [46].
Therefore reflection would be mainly produced on an almost clean surface in spite of the large amount of incoming O atoms reaching the surface. Our previous microkinetic model confirmed as well a very low coverage at quasithermal conditions [5]. However, the present simulations are performed over an O-precovered surface (coverage close to 100% but only in the initially sampled (1 ´ 1) unit cell, allowing the atomic adsorption only in the empty neighboring cells), which shows a different behavior due to the influence of adatom. This coverage effect was partially observed in previous theoretical simulations of O impinging several models of graphene/graphite surfaces [2,13]. According to this, the coverage of HOPG samples during the experiment [2] should be small and the samples are expected to remain essentially clean after experimental O/O2 mixture exposure. Other factors that could also affect the O scattering distribution would be the inclusion of the breaking of C sheet (i.e., vacants favor O adsorption) produced by the formation of CO and CO2 and its deformation during the total process, not included in the present study.
The O2 molecules formed via ER reaction have also been analysed showing an angular distribution ( Fig. 10b) peaked at lower angles (peak at ~ 40°) than reflected O atoms, close to the specular angle.
There is as well some experimental data [2]  In order to make a more reliable comparison with experimental data [2], an extensive QCT study of O2 molecules impinging a clean surface is also in progress [68], which tends to produce angular distributions peaked at higher angles than those obtained for ER distributions, with a good agreement with the experiments.
The previous molecular dynamics simulations for hyperthermal collisions of atomic oxygen over some graphene models [2,13] confirm that O2 is formed via a direct ER mechanism and that this process along with epoxide formation and inelastic atomic scattering are the main channels. These simulations do not control the surface temperature (e.g., this can increase up to 4,500K [13]) and only run 100 trajectories (i.e., large statistical errors) for normal incidence and Ecol = 5 eV. Thus, a direct comparison with the present results is not accurate. However, a qualitatively comparison can be intended with respect to their models of graphene sheets functionalized with epoxides, which are closer to our model of unit cell, although their atomic coverages were also increased during the simulations.
In general, it is observed that ER process is the most important one (e.g., 72 % for model IV, 30 % for model V [13] or 76-79 % for model IV [2]) to be compared with 43,8%, here calculated at the simulated experimental conditions. They also observed that O reflection was only observed for larger models (e.g., 17 % for model V [13]) where sheet deformation was important. Our calculated value is around 31.9 %. The QCT larger values (e.g., comparing with model V) are consistent with the sampling of a (1 ´ 1) unit cell with an O-preadsorbed, which prevents a large adsorption, hence increasing the other channels, as ER or atomic reflection.

V. Summary and conclusions
A new PES for the interaction of atomic oxygen with an O-preadsorbed graphite (0001) surface has been constructed. The analytical FPLEPS PES is based on a set of 563 DFT data points calculated over high-symmetry sites centred in the ER channel. The symmetry of the graphite surface is fulfilled for the two-atom interaction by using appropriate periodic functions to describe the angular dependence and the X-Y dependence of the potential. The agreement between FPLEPS and DFT data is quite good, with a RMSD of 0.120 eV for energies lower than 1 eV over reactants.
The QCT method has been applied for the dynamics study of O collisions over O-preadsorbed graphite (0001) surface as a function of collision energy (0.01 -2.0 eV), surface temperature (100-1,500 K) and incident angle (0°, 45°). Moreover, quasithermal and hyperthermal conditions were also simulated and compared with available experimental data.
For normal incidence, ER reaction and O reflection are the main processes in agreement with earlier experimental and theoretical studies, becoming the ER more important than the O reflection for collision energies higher than 0.2 eV. For off normal incidence, ER recombination is always lower than O reflection. All other processes are almost negligible (e.g., atomic desorption) for both normal and off normal incidence.
Polar scattering angle distributions of atomic oxygen are shifted to higher angles compared to collisions with a clean surface for normal incidence, while for off normal incidence this shift corresponds to smaller angles, showing that preadsorbed oxygen affects these angular distributions.
Molecular oxygen formed via ER reaction is translationally and internally excited and presents polar scattering angle distributions centred at low angles, which are rather similar for all conditions studied.
The increase of surface temperature introduces small changes in almost all studied properties.
Moreover, the study for reactants at quasithermal conditions is in agreement with the main features observed for studies at low initial collision energies.
Polar scattering angle distribution for O colliding with an O-preadsorbed surface at hyperthermal experimental conditions at θv = 45°, Tsurf = 503 K shows a peak at low angles (~35°) far from experiment data (~62°), but the calculated polar scattering angle distribution of atomic oxygen colliding onto a clean graphite surface matches well the experimental one. [46] However, the comparison between theoretical and experimental distributions cannot be straightforwardly done because the experimental HOPG surface coverage is unknown, and we show that coverage effect is significant. According to our results, we could conclude that the coverage for the HOPG samples during the experiment should be small, so they would remain essentially clean after the experimental conditions reported.
As possible improvements to the present study we would mention the possible inclusion of the breaking of C layer to produce CO, CO2 and vacants and also the layer deformation during the total process.
distributions of atomic and molecular oxygen in the experiments for hyperthermal O( 3 P)/O2 colliding over HOPG. We are also grateful to Dr. Fabio Busnengo for sending us the GLO subroutines and for useful discussions about this method.    ) for B1, B2, B3, T1, T2, H,   T1H1, T1H2, T1H3, T2H1, T2H2, T2H3, respectively. The minimum triangle shaded was used for the atomic interaction with the clean surface [46]. Two layers of graphite are depicted: black circles (1st layer) and gray circles (2nd layer).