Classical dynamics study of atomic oxygen over graphite (0001) with new interpolated and analytical potential energy surfaces

Two adiabatic potential energy surfaces (PES) based on density functional theory data are constructed for the interaction of atomic oxygen with graphite (0001) surface: an analytical FPLEPS PES and an interpolated Modified Shepard one. A classical trajectory study has been performed for the two PES for different initial conditions: collision energy (0.1 £ E col £ 1.3 eV), surface temperature (100 £ T surf £ 900 K) and two incident angles ( q v = 0 ° , 45 ° ), and also for thermal conditions (T = T Oxygen = T surf = 300-1,500 K). In addition, hyperthermal experimental conditions corresponding to a hot atom distribution (<E col > = 5. 2 eV) were also considered. All the properties studied for the two PES were in close agreement in almost the major part of the explored conditions, although some differences were obtained for low E col due to the presence of a physisorption minimum in the MS PES that was not included into the FPLEPS one. The adsorption process occurs mainly over bridge sites. Adsorption probabilities are lower than reflection ones in practically all the conditions explored and increase quickly with E col until a maximum and then decrease smoothly. Polar is mainly from the atom to the surface and increases when initial collision energy does. experimental velocity and translational distributions of O/O 2 beam species, used for the experiments with hyperthermal O colliding over HOPG.


Introduction
Carbon-based materials have received considerable attention for their use in spacecraft applications. In particular, organic thin films, polymers and carbon fibre reinforced composites among others are used as a thermal protection system (TPS) in some part of Space Shuttle due to their light weight and high strength. In low Earth orbits (LEO), at approximately 160-2,000 km above Earth's surface, atomic oxygen is one of the main present species and can collide with the spacecraft surface at very high translational energies (i.e., about 4.5±1.0 eV due to the orbital velocity of around 7.5 km/s) [1]. Some studies obtain information about the damage (i.e., etching, degradation,...) that produces these highly energetic collisions over the materials in order to improve them. In addition, a close understanding of the reactions that occur between graphite and hyperthermal O atoms is also of great interest to the rocket plume research community.
Experiments are performed with beams of O atoms (although some beams are formed by significant concentration of both O and O2 in their ground state) that impinge over highly ordered/oriented pyrolitic graphite (HOPG). Some of the experiments [2,3,4] were focused on the study of the chemical reactivity and morphological evolution of HOPG upon exposure to 5 eV oxygen atoms (hyperthermal conditions). After exposure to atomic oxygen, the graphite surface becomes functionalized with epoxide groups that can migrate across the surface. Moreover, molecular oxygen is expected to scatter inelastically from the surface since it is only able to physisorb over a pristine surface by a small barrier. Neither 5 eV oxygen atoms nor 10 eV oxygen molecules are expected to cause sputtering from a pristine unfunctionalized surface [5] at least if no defect sites are formed.
Some theoretical works about the interaction of atomic and molecular oxygen with graphite [6,7,8,9,10,11], graphene [12,13] and carbon nanotubes [8,14,15] were mainly concerned with the adsorption over several types of surfaces for different coverages. Atomic oxygen becomes chemically adsorbed on a bridge position between two adjacent C atoms of graphite or graphene, forming an epoxi group with an adsorption energy in the range of 0.95 -3.2 eV [6,7,8,9,10] 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 energies increase considerably (4.5 -5.7 eV) [10]. Another important feature of the atomic oxygen adsorption over graphite is the easy diffusion observed (energy barriers of 0.38 eV [6] and 0.36 eV [7] for a coverage of 12.5 % and 3.12 %, respectively), which are lower than the obtained for a graphene sheet (0.58 eV for a coverage of 16.7 % [12]), that would imply a noticeable mobility of isolated adsorbed atom. In addition, we recently proposed a microkinetic model for O/O2 over graphite [6].
This model used rate constants obtained from DFT and standard transition state theory for several processes occurring when O/O2 mixtures collide over graphite surface showing a very low steadystate atomic coverage (q < 0.5%) and also very low atomic recombination coefficients at thermal conditions (300-1,000 K).
Two recent works have simulated hyperthermal normal collisions of O over graphene [2,16] and graphite [16] by using DFT (PBE/DZP) data and the ReaxFF reactive force field, respectively, but without inclusion of a thermal bath. The surfaces were functionalized with preadsorbed O atoms or single vacancy defects. One of these dynamical simulations allows very long time integrations (e.g., 150 ps [16]), being the consecutive O collisions with the previous surface produced every picosecond, which increases also the system size by one atom per collision. In these conditions, it was observed that O adsorption (i.e., epoxide and carbonil formation) and O2 formation via an Eley-Rideal (ER) mechanism were the main processes (t < 45 ps) but inelastic collisions were very low.
At longer times CO/CO2 formation could also be observed and C sheet was both broken and deformed.
The main goal of this work is to provide new potential energy surfaces (PES) for studying the O collisions with a pristine graphitic (0001) surface and to obtain dynamical information of this system by means of classical trajectories (CT), which can also be compared with some experimental data [2]. Several PES are in the literature for studying the hydrogen recombination over graphite and graphene surfaces (i.e., H(g)+ H-graphite [17,18], H(g) + H-graphene [19]).However, up to our knowledge, there are not available analytical or interpolated PES for describing the atomic interaction of oxygen with a clean graphite surface, unless an earlier PES derived assuming O/ethylene data [20] and used for calculating the erosion yield (i.e., loss function) of graphite by O atoms.
This paper is organized as follows: Section 2 provides a brief detailed description about DFT and CT calculations; Section 3 describes the PES construction and their 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.

Theoretical methods
In this section, and for brevity purpose, only the main details of density functional theory (DFT) and classical trajectory (CT) calculations are given. The description of analytical and interpolated methods used in the potential energy surface (PES) construction will be done in Sec. 3.

DFT calculations
Periodic DFT calculations have been performed using the Vienna ab initio simulation package (VASP) [21,22,23,24], that uses plane waves basis set. A detailed description of the DFT calculations was done in a previous work [6], thus only a few details are given here. The calculations are mainly based on the generalized gradient approximation (GGA) with the revised Perdew-Burke-Ernzerhof functional (RPBE) [25]. The projector-augmented wave (PAW) technique within the frozen core approximation has been used to describe the electron-core interaction [26,27]. An energy cut-off of 550 eV has been used in the plane-wave expansion after several tests to confirm the energy convergence of the calculations. Integration over the Brillouin zone was performed by using an (11 ´ 11 ´ 1) k-point mesh by means of the Monkhorst-Pack method [28] for slab calculations using a (2 ´ 2) supercell. The energy convergence in the electronic selfconsistent procedure was maintained below 10 -6 eV for all geometrical calculations performed under the rigid slab model. The vacuum used in the slab model was around 15 Å, which was large enough to prevent significant interactions between periodic images. Due to the important role that spin plays in atomic oxygen (i.e., triplet in its ground state), all calculations were spin-polarized [6,8,10].
Graphite is formed by a layered structure of six-member carbon rings, stacked in an ABAB sequence with a weak interlayer binding compared to the much stronger binding within the layers, (Fig. 1a). The experimental [29] lattice parameters for graphite structure are a = b = 2.456 Å and c = 6.696 Å. The theoretical values obtained in our previous work [6] reproduce the a and b parameters with a value of 2.468 Å but overestimate the c one. Thus, we decided to use the experimental c lattice parameter throughout the DFT calculations since this assumption led to a good description of the rest of bulk properties. Three different sites were considered in the case of atomic oxygen adsorption: on top of a surface C atom, which is either above another C atom (T1) or above an hexagon hollow (T2) of the second layer, and over a bridge (B) between the two nearestneighbour carbon atoms (Fig. 1b). The two top sites (T1 and T2) are practically identical since the spatial difference between them is in the layer below, which is located at a distance of 3.348 Å, too far to produce significant differences in the interactions. We checked also that the hollow site (H) of the hexagonal ring of the first layer did not present any significant adsorption in agreement with previous studies [30,10]. The adsorption energies obtained for a (2 ´ 2) supercell show that the most stable location for O-adsorption over a pristine surface is onto the B site (Eads = 0.660 eV) while for T1 and T2 sites the adsorption energy is practically the half (Eads = 0.308 and 0.281 eV, respectively). Moreover, we have observed that spin-polarized calculations with free magnetization were almost coincident with the lowest spin state energies in each region; at large distances to the graphite (e.g., d(O-surface) > 2.0 Å) the lowest energy corresponds to the triplet state but at shorter distances the singlet state has a lower energy. Consequently, we assume that the evolution of the system is electronically adiabatic, which requires at least one transition (i.e., for adsorption) between states of different spin (intersystem crossing, ISC). In spite of the possible importance of nonadiabatic effects in some gas-surface processes (e.g., for metals, [31]), in the present system the different earlier experimental and theoretical works can support the adiabatic approach. Thus, an ab initio direct Ehrenfest dynamics study with time dependent DFT (TDDFT) for collisions of O( 3 P) with graphite clusters confirms the change from triplet to singlet character at the top and bridge sites at the closest approach during the collision [32]. This triplet/singlet ISC was also assumed to be important in a dynamical study of graphite erosion by atomic oxygen [33]. A recent molecular dynamics study of O/graphene using ReaxFF reactive force field [16] shows that epoxide formation is the main process, which is also experimentally observed for O collisions over graphitic surfaces [34,35] producing graphite and graphene oxides.
DFT curves calculated over a specific site for several Z distances (

CT calculations 8
A dynamical study of the adsorption and reflection for the atomic oxygen over graphite was done by means of the classical trajectory method [37] (CT). Several initial atomic collision energies between 0.05 eV £ Ecol £ 3.0 eV were considered for both FPLEPS and MS surfaces at two incident angles, θv = 0° (normal incidence) and 45°, of the initial velocity vector respect the normal vector to the surface. All these conditions were explored assuming rigid slab model, without the effect of temperature and also including a Generalized Langevin Oscillator (GLO) [37,38,39,40] model to simulate the surface temperatures (Tsurf) in the range of 100 K £ Tsurf £ 1,300 K. In addition, thermal initial conditions for both O atom and surface were investigated within temperature range between 300 K £ T £ 1,500 K for initial θv = 0° and 45° approaching angles (fv angle was uniformly sampled between 0 -360°). Finally, a hyperthermal atom distribution with graphite surface at Tsurf = 503 K was also studied in order to compare with experimental data available [2].
Initial incoming oxygen position (X, Y) is randomly selected along the (1 ´ 1) unit cell meanwhile initial Z position is set to 6.5 Å where the interaction with surface is negligible.
Classical trajectories were calculated with qctsurf program developed in our group, which integrates the Hamilton equations of the system (O + GLO) using the Beeman algorithm. The time step used for both surfaces was 5´10 -17 s, which ensure 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 1.5´10 -12 s, long enough to ensure proper interaction of the oxygen atom with the surface.
In order to classify the trajectories in one of the two possible exit channels (adsorption or reflection) the following criteria were taken into account. For a trajectory to be considered as reflection, the Z coordinate should be higher than 7.1 Å (i.e., no atom-surface interaction) and the direction of the velocity vector, should be pointing to the vacuum. On the other hand, for a trajectory to be considered as adsorption, the Z of the incoming atom has been considered to be lower than 2.2 Å (a value close to the one reported for the adsorption barriers) [6] and the number of total rebounds with the surface should be at least 8, irrespectively of the final position. The minimum number of total trajectories (NT) calculated for each condition was around NT = 10,000 which, for example, leads to a standard error of 1 -2% in adsorption probability.
To account for energy exchange and dissipation a GLO model [37][38][39][40] was used, including the surface and ghost particle movement into the Hamilton equations. In the present work, as initial oscillator frequencies (wi,x, i= 1,2,3) were used the values derived from bulk (TD) or surface Debye temperature (Ts,D) of graphite (i.e., Ts,D= 350 K [41] and TD= 423.6 K [42], 393.9 K [43]). Thus, values close to 10 -3 for wD and to 5´10 -4 for wi (in atomic units of time -1 ) were calculated. These values were used as guess frequencies to fit the experimental mean square of perpendicular displacements of graphite surface atoms as a function of surface temperature [44].
The agreement with experimental data is quite good for surface temperatures between 150 -300 K. The optimal parameters are wi,x = wi,y =10 -3 and wi,z =3.4 ´10 -4 a.u. for i=1,2,3. It was checked that changes on the effective mass [e.g., from 12 (one carbon atom) to 60 amu (5 carbon atoms of the unit cell)] produce negligible effects on the quality of the fit of the experimental values. In fact, the incoming atom can interact with more than one surface atom (Armand effect [45]) during the collision due to the finite size of the projectile, which can justify even more the use of an effective mass. The value of 60 amu along with the chosen friction constants (gg,x = gg,y = gg,z = 4.0 ´10 -4 a.u.) were good enough to reproduce the experimental hyperthermal data, shown in section 4.
Within the framework of the FPLEPS, the 3D PES expression for atomic oxygen over graphite can be developed by the equation, where X, Y, Z are the coordinates of the oxygen atom. The first term is a modified Morse potential, where D and a are the Morse parameters and Z eq is the equilibrium Z distance of adsorbed oxygen. a parameter is expressed by the function, where ah and al functions ensure a good description of both the attractive and the repulsive part of the atom-surface potential, which are smoothly connected using a switching function (f), According to this, a = al for values of Z < 1.20 Å (the repulsive branch) and a = ah for values of Z > 1.66 Å (the attractive branch).
The second term introduces more flexibility to the Morse potential to be able to reproduce several adsorption energy barriers at different sites of the solid surface. In earlier works (e.g., N2/W [49,50,51]) this term was not necessary. Thus, the used J(X,Y) and c(X,Y) functions are, The second term in equation (1) varies from 0 to 2l when Z changes from zl to zh.

11
There are a total of nine parameters (D, Z eq , a0, a1, a0 ', a1 ', zl, zh, l), whose (X,Y) dependence is introduced via a Fourier expansion representing the symmetry of the crystal surface.
For each Pj parameter (e.g., P1 º D, P2 º Z eq ,..) a Fourier series up to the 5 th order is used for graphite (0001) surface, (8) where Cji (i = 1 to 5 and j = 1 to 9) 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) Six selected high-symmetry sites (Fig. 1b), namely T1, T2, B, H, T1H and T2H, inside the minimum triangle that generates the entire unit cell, were considered for DFT calculations. Thus, a total of sixteen Z values for each curve were calculated in the range of 0.2 £ Z £ 3.1 Å (i.e., a total of 96 DFT points). Then, a least squares fitting procedure was performed independently for each of the 6 DFT curves, which correspond to the high-symmetry sites. The optimal parameters are shown in Table 1 and they reproduce quite well the calculated DFT curves ( Fig. 2a and 2b). The Cji Fourier coefficients, which will guarantee the (X,Y) dependence in other sites of the surface, were obtained by solving a system of linear equations, and are given in Table 2.
Analytical derivatives were also obtained for this PES in order to get a more accurate and faster integration of the Hamilton's equations.
In spite of the LEPS, PLEPS or FPLEPS names are frequently applied to 3-atom systems (gas phase) or 2-atom/surface systems, we have maintained the name FPLEPS for the O/graphite PES because this is a particular case from the most general FPLEPS surface that we are also constructing for the O2/graphite system.

Interpolated MS PES
We have also constructed an interpolated PES for O/graphite (0001) system using a modified Shepard scheme, which was initially developed by Collins and co-workers for gas-phase potential energy surfaces [52,53,54,55,56] using internal coordinates. Lately, this methodology was applied also to some gas-surface reactions (H2/Pt(111) [57,58], H2/Pd(111) [59]) with some difficulties (e.g., derivative discontinuities in the dynamics). We have developed a PES based on Cartesian coordinates. The energy is evaluated as a weighted sum of Taylor series expansions centered about a number of calculated DFT points (Ndata), (10) where Ti (X,Y,Z) is a second order Taylor expansion of the potential centered on the i-th data point being (X,Y,Z) º (q1,q2,q3). The weight function (ni) is calculated from an inverse distance power between (X,Y,Z) and (Zi,Yi,Zi) points, and then is normalized (wi) by summing over all Ndata points,

13
The selection of p parameter was tested for different values within the range 4 ≤ p ≤ 12, checking also the energy conservation during some trajectories, obtaining an optimal value of p = 8. Using the same checking procedure, a threshold or tolerance for the normalized weight was stated at wtol = 10 −8 . Thus, configurations with a smaller weight (i.e., wi ≤ wtol) were not considered in the total sum of equation (10) to save computer time without reducing the energy accuracy (Neffective < Ndata).
In order to describe properly the asymptotic regions for Z > Z1= 3.7 Å, the following equation was used for energy calculation, (13) where V ass is the asymptotic energy (e.g., O + slab) and f(Z) has the expression, This latter function ensures a soft connection for energy and derivatives in these regions.
In spite of MS interpolation does not need a regular grid of DFT data points, we have preferred to use an initial regular (u,v) grid of DFT points (Ndata = 399), calculated inside the minimum triangle area (green zone in Fig. 1b) for 0.9 Å £ Z £ 3.7 Å. Equivalent configurations were obtained by using symmetry operations for describing not only all the (1 ´ 1) unit cell (i.e., mirror planes) but also half of the surrounding cells (i.e., lattice translations). Thus, a total of 9,156 configurations were used, of which 2,289 were inside of the initial unit cell. DFT energies, analytical gradients and numerical second derivatives were computed for all data points. The initial PES data set was iteratively grown by adding 66 extra DFT points for three initial collision energies (Ecol = 0.3, 1.0 and 3.0 eV), which gave place to a final value of Ndata = 465. Two criteria were alternatively applied in each cycle to select the new points [57]. The first one (h-weight criterion) chooses the configurations most regularly visited by the trajectories (i.e., configurations with the largest h values) and described with a few number of DFT data points. The second one (variance sampling criterion) selects points in the regions where the calculation of the energy seems to be less accurate (i.e., configurations with the largest variances).
Analytical derivatives were also obtained for the MS PES as in the previous FPLEPS.

Analysis and comparison of both PES
A first analysis of the FPLEPS surface can be performed looking at the accuracy of the fitted DFT grid of points. Thus, the calculated root mean square deviation (RMSD) respect to the 64 DFT energies used for the fit, corresponding to 1.1 Å £ Z £ 3.1 Å (the most repulsive were not included), was of 0.034 eV, which is small enough when it is compared with expected DFT uncertainties.
Moreover, we have also calculated the RMSD for the DFT grid used latter in MS PES, excluding both the most repulsive points and those included in the same fit. Thus, a value of 0.095 eV for the 328 DFT points is achieved, which is quite good because these points were not used for the fit. and H), respectively. They compare reasonably well. A similar behaviour is found in Fig. 2c for the MS curves. Nevertheless, FPLEPS curves are much more smoother than the interpolated ones.
Moreover, one additional curve corresponding to the hollow of the centre of the triangle in Fig. 1b (hereafter H2 site), where there are not DFT points, is also calculated for both PES, showing a very similar shape. Table 3 summarizes the minimum and the energy barrier properties at five sites for DFT, FPLEPS and MS curves. It is important to note that the absolute energy used as the asymptote for the FPLEPS surface, differs slightly from the one used in MS PES. In all cases, the zero of energies was obtained for a supermolecule calculation considering the oxygen atom far away from the surface although the absolute value for the MS PES was 0.052 eV higher. As above mentioned, the MS curves (not the FPLEPS ones) tend to show a somewhat vdW minimum around 2.75 Å with adsorption energy lower than 0.06 eV, which is more stabilized once dispersion energies are included. Nevertheless, this physisorption minimum is only important at very low collision energies, whereas the energies implied in the applicability of the TPS are by far much higher. MS PES is built by using a larger number of DFT geometries, so it is more appropriate for describing the vdW minimum zone. The FPLEPS could describe also this minimum by introducing additional terms (i.e., Gaussian functions) in equation (1) that would give more flexibility in the entrance channel.
The description of B, T1 and T2 sites is similar for FPLEPS and MS surfaces although the B adsorption minimum is slightly more stable in MS PES. This is not only due to the reference asymptote difference but also to the worse quality of the MS first and second derivatives around this minimum. It is shown that adsorption is an activated process as there are energy barriers in all studied sites. The presence of this activation barrier arises from the ISC between the triplet and singlet states in the same manner that MRCI method provided an activation barrier of 0.73 eV for the interaction of oxygen with C2H4 [60]. Table 3 also shows the site with the lowest adsorption barrier (hereafter TS site). This is between T1 and T1H sites or close to B site for FPLEPS and MS PES, respectively, although in both cases their energy barriers are only somewhat lower than for B site. The similarity of T1 and T2 curves ( Fig. 2) and data (Table 3) is consistent with the values found in previous DFT calculations [6].
Several energy contour plots at two fixed oxygen Z distances to the surface as a function of X and Y values are presented in Fig. 3 in order to facilitate a complete view and comparison between both surfaces. In general, good agreement between both representations is found. Clearly, the bridge sites are the most attractive regions for atomic adsorption and the hexagonal hollow sites the most repulsive ones. Therefore, it will be expected that O tends to be adsorbed close to B sites and subsurface penetration through hollow sites will be forbidden.
MS PES shows wiggly equipotential lines due to the effect of low accuracy on the second derivatives.

Probabilities
Figs. 4a and 4b show that in absence of a thermal bath, oxygen reflection is more important than atomic adsorption at all energies and incident angles (e.g., θv = 0° and 45°) explored for both PES (FPLEPS and MS). It is worth noting that the main difference between PES is found at low energies (i.e., Ecol = 0.05 eV), especially for off normal incidence where the adsorption probability is even higher than the one obtained for the reflection process in the MS surface. This effect is due to the presence of the small physisorption minimum in the MS surface, not reproduced in FPLEPS one.
Moreover, there is threshold energy around 0.2 eV in the adsorption probability with the FPLEPS PES for normal and off normal (in this case the threshold corresponds to perpendicular energy) incidence. These threshold values are in agreement with the lowest energy barrier found in the PES (Table 3), although the energy barrier heights depend on the position of the incoming atom. Thus, when the atom overcomes the barrier to the potential well (Z< 2.0 Å), a part of the perpendicular kinetic energy can be transformed into the parallel one, consequently favouring the trapping of the atom. For a small increment of the collision energy the chemisorption probability can increase until reaching a maximum value at around 0.325 eV (i.e., 27% of adsorption probability for both PES), because at these higher energies the incoming atom is able to overcome not only the minimum adsorption barrier but also higher ones located in close surface regions (Table 3). However, at much higher collision energies the trapping finally decreases as could be expected.
For off normal incidence (Fig. 4b), it is observed a clear increase of the reflection at all perpendicular energies with a maximum on the adsorption probability around 0.4 eV (1.9%) and 0.5 eV (2.5%) for FPLEPS and MS surfaces, respectively. These adsorption probabilities are much lower than for normal incidence and show also the same threshold energies.
Figs. 4c-4f present the reflection and adsorption probabilities once introduced the surface temperature (Tsurf = 100 K and 900 K). Adsorption probabilities show now higher values than without thermal bath at the same energies for both incident angles and both PES. In the case of normal incidence, adsorption is the main process in the collision energy range 0.35 -0.95 eV and 0.50 -0.90 eV for FPLEPS and MS PES, respectively, overcoming the reflection probability, which was not observed without thermal bath. The thermal bath allows a noteworthy decreasing in the perpendicular energy due to an energy transfer to the graphite surface explaining the high increase of the adsorption probability. For off normal incidence, the transfer from perpendicular to parallel is limited because we have initially parallel energy (like a communicating vessels system) and the adsorption probability is lower. Adsorption threshold energy is observed to occur at the same initial perpendicular collision energies when the thermal bath is added, because this latter is not only giving extra stabilization by exchanging energy with oxygen atoms but also it is allowing that atoms with smaller collision energy can overcome the energy barrier. The effect of energy exchange between atomic oxygen and the graphite surface is more important as the surface temperature is increased but mainly for lower collision energies. At higher collision energies the small increase of surface temperature from 100 to 900 K limits the atom-surface energy exchange, as a result decreasing a little the adsorption probability.
The MS adsorption probabilities show a minimum around 0.2 eV not found in FPLEPS ones, with a noteworthy adsorption at lower collision energies. There is a predominant chemisorption show the final positions of the adsorbed atoms for normal and off normal incidence and Ecol = 0.5 eV, which are mainly located over the bridge sites and their surroundings areas, whereas no adsorption is found over the repulsive areas (e.g., hexagonal hollow sites) that correspond mainly to the initial positions of the reflection processes, as predicted by the previous theoretical calculations [6]. It is worth noting that trajectories initially started inside the central (1 ´ 1) unit cell can also become adsorbed over bridge sites of adjacent cells. This could be explained considering that oxygen atoms can easily diffuse over the graphite surface before becoming adsorbed, in agreement with some experimental [3,4] results and also with the low diffusion barriers reported in previous DFT studies [6].  Table 3), showing that reflection happens for distances similar or larger than those of the adsorption barriers, leading to close specular scattering angles. For higher collision energies the incoming atom is able to get closer to the surface, overcoming easily these barriers and achieving more repulsive regions of the PES, which can justify the separation from the specular angles.

Polar scattering angle distributions
The surface temperature effect is small between 0-900 K; the increase of Tsurf broadens slightly the peaks. Moreover, we have carried out a deeper analysis of the scattering distributions at several collision energies (Ecol = 0.25, 0.5, 1.3 eV) for Tsurf =300 K and θV = 45º looking also the parallel, perpendicular and total energy transfers from the incoming atom to the surface along with the minimum approach distance to the surface (i.e., <Zlow>) around the maximum (± 5°) of the scattering distributions (I(θV')). For an approaching angle of θV = 45º, the initial perpendicular energy is always half of Ecol. Several trends are observed regarding both the shift and shape of 20 I(θV') distributions. It is worth noting that larger θV' angles than the initial ones can be originated by a decrease of the perpendicular energy or/and an increase of the parallel energy. Thus, from 0.25 to 0.5 eV the peak shift (from 47° to 59°) is only produced by a perpendicular to parallel energy transfer with a negligible contribution of the solid surface. Since the <Zlow> values are larger than 2 Å, the perpendicular energy is still lower than the barriers (Table 1, Fig. 2). At 1.3 eV the distribution is much wider and peaked at 35°. This change is mainly produced because now the parallel energy is vastly transferred to the solid surface, achieving shorter distances (<Zlow> = 1.1 Å).
A few differences are observed between FPLEPS and MS PES distributions. There is a general tendency of MS PES to produce polar scattering angle distributions shifted to large angles, possibly originated by the long range (i.e., physisorption well) PES differences. Fig. 7 shows the distributions of the energy transferred by the reflected atoms to the surface for normal incidence, several initial collision energies (Ecol = 0.25, 0.50 and 1.30 eV) and two surface temperatures (Tsurf = 300, 900 K) for FPLEPS and MS PES. This energy is calculated from the difference between the final and initial atomic energies (i.e., ∆E = E'col− E col), being negative if the atom transfers energy to the surface. In general ∆E < 0 and becomes more negative as the initial collision energy increases. Nevertheless, it can be observed for a small number of trajectories a positive energy exchange, which means that in some cases, the surface is giving energy to the impinging atom. This effect that increases with temperature is mainly observed at low initial collision energies.

Energy distributions
Some differences are observed in these distributions between FPLEPS and MS PES at higher surface temperatures and collision energies (Fig. 7b). MS peaks are slightly shifted to more negative ∆E values and become narrower. Therefore, it seems that MS PES favours more in such conditions the exchange of energy from the incoming atoms to the surface. This fact could help to understand the differences in Fig. 6 for both PES. If there is more energy exchange (MS PES) the atoms interact more strongly and can lose much more the memory of the initial conditions, leading to broader angular distributions. Moreover, the higher adsorption probability observed with MS PES (Fig. 4) at high collision energies can be explained also according to this more important energy exchange.

Thermal and hyperthermal initial conditions
Thermal probabilities at several temperatures (i.e., T = TOxygen = Tsurf) and at two incident angles (θv = 0° and 45°) for both PES were also computed to facilitate also the possible comparison of CT rate constants for adsorption process with previous TST ones [6] and with future experimental measurements. Results are plotted in Fig. 8, which shows very low adsorption probabilities that increase as temperature does in all the range explored (300 £ T £ 1,500 K) for FPLEPS PES.
Normal incidence adsorption probability (Fig. 8a) is slightly higher than the one for off normal incidence (Fig. 8b). The general trend of these probabilities is close to the one found for low initial collision energies, including surface temperature (Figs. 4c-4f) and also is comparable to the previous results of our microkinetic modelling [6], which predicts a very low atomic adsorption at thermal conditions (300-1,000 K).
FPLEPS and MS probabilities are almost coincident at high temperatures (i.e., T > 1,300 K). At lower temperatures the MS adsorption probability decreases when T augments, in a similar way as for the increase of collision energy until 0.2 eV, due to the physisorption minimum effect. When increasing the temperature, the physisorption process can decrease mainly by two reasons: firstly, because the initial collision energy is increasing, so atoms have too much energy to be stabilized on the physisorption well, and secondly, because a higher surface temperature allows a higher energetic exchange from the surface to the oxygen atom, which can give energy enough to the atom to go back to the gas phase. On the other hand, when temperature is high enough, some of the atoms will have energy sufficient to surmount the adsorption barrier and to lead to chemisorption. Thus, the increase in adsorption probability from 1,300 K to 1,500 K will be produced by the increment of the chemisorption.
The polar scattering angles observed at several temperatures are comparable to the scattering angle distributions observed at low initial collision energies (Figs. 6c-6f).
Polar scattering angle distribution was also computed at the particular experimental hyperthermal conditions recently measured [2]. Thus, 450,000 trajectories were computed for both PES at θv = 45°, Tsurf = 503 K, sampling the same experimental initial O atom distribution, with an average value (<Ecol>) of 5.2 eV. Fig. 9a compares the calculated angular distributions with the experimental one, which was measured for scattered atoms in the same plane as the initial plane formed by the normal to the surface and the O beam line (i.e., in-plane measurements). Scattered atoms are considered in-plane when the final plane was within 1° (angular tolerance) respect the initial one. Calculated in-plane distributions for MS and FPLEPS PES are similar and close to the experimental curve. A smaller angular tolerance tends to narrow the calculated distributions. On the other hand, out-of-plane contribution produces broader distributions (Fig. 9b) although peaked at similar angle. In addition, the number of scattered atoms out-of-plane is larger by far than the inplane ones (ca., 98.7% out-of-plane and 1.3% in-plane for FPLEPS and MS with an angular tolerance of 1°). Therefore, experimental measurements out-of-plane should report even a broader angular distribution.
The distribution obtained in Fig. 9a is super-specular, while for lower energies (Ecol = 0.25, 0.50, 1.30 eV, Fig. 6) the distributions were peaked at specular or even lower angles. It is important to note that when Ecol is increased sufficiently over the barrier regions (e.g., 2.0, 3.0, 5.2 eV) we observe that there is a strong perpendicular energy transfer to the solid surface together with a perpendicular to parallel energy transfer (more evident when Ecol is increased and clearly at 5.2 eV), which produces peaks at larger angles, changing from 55° to 65°. Summarizing, the different regions (i.e., corrugation) of the PES available at each collision energy are clearly responsible of these changes in the I(θV') distributions.
These peaks are very close (θv' = 65° for the theoretical results and θv' = 62° for the experimental one), although experimental distribution is somewhat narrower than the theoretical ones as above mentioned. However, some consideration should be made. Secondly, the breaking of the C sheet (produced by the formation of CO and CO2) and its deformation [16] could also alter the experimental scattering distribution. Finally, the use of a better thermal bath with a complete slab movement (full PES) could be necessary to simulate better these experimental high energetic conditions, which would also provide somewhat larger adsorption probabilities than the ones we obtain for both PES in these conditions (ca. 0.5%). Previous dynamical simulations [16] confirm that the inclusion of the second C layer (graphite model) acts a heat sink and also delays the creation of defects respect the graphene model. In fact, the present CT simulations (one O atom collision) are not directly comparable with these ones [16] (ca. 150 O atoms colliding per trajectory) because in our study the surface is always clean and pristine whereas in this another work the surface becomes altered during the much longer simulations (ca. 1 ps vs.

24
150 ps) even when the simulations start with a model of pristine either graphene or graphite surface.
The distributions obtained by using both PES are in good agreement as could be expected because at these extreme collision conditions the small PES differences will be insignificant.

Summary and conclusions
Two different potential energy surfaces for the interaction of atomic oxygen over graphite (0001) surface have been constructed. The analytical FPLEPS PES is based on a reduced set of DFT data over high-symmetry sites, and the interpolated Modified Shepard one (MS PES) uses a much larger set of DFT data, that were calculated in the present work for this goal.
Classical trajectory method has been applied to study the adsorption and reflection processes as a function of the collision energy (E col = 0.1-1.3 eV), the surface temperature (T surf = 100-900 K) and the incident angle (θv = 0°, 45°). Moreover, thermal and hyperthermal conditions were also considered to compare with available experimental data.
In general, the results obtained by the two PES used are in close agreement in practically the major part of the explored conditions. Nevertheless, some differences were obtained for low collision energies due to the presence of a physisorption minimum in the MS PES that was not reproduced in the FPLEPS one. The description of the system by means of the interpolated PES could be considered more accurate than the analytical PES since a large amount of DFT information was included (not only energies but also first and second derivatives) in order to describe the system. It is worth noting that MS surface matches completely the DFT information whereas the fitted FPLEPS has an associated root mean square deviation.
Adsorption probabilities are lower than reflection ones for almost all initial conditions studied.
Adsorption increases quickly with collision energy until a maximum and then decreases smoothly.
This behaviour is coherent with the similar adsorption energy barriers found on several sites of the 25 graphite surface. Especially, the bridge sites correspond to the most favourable minimum energy path for the adsorption. Only at collision energies lower than 0.2 eV there is a clear difference between both PES. FPLEPS shows no adsorption whereas MS PES shows a significant adsorption that increases at even lower collision energies, due to the physisorption contribution. The use of several incident angles and surface temperatures produces only small changes in these probabilities. Polar scattering angle distributions are similar for both PES, although with some minor discrepancies at low collision energies. The peak of these distributions for normal or off normal incidence is mainly observed around their specular position, although the increase of the collision energy broadens these distributions and moves them at larger angles. The increase of surface temperature enlarges slightly the peaks.
The atom-surface collisions mainly transfer energy from the atom to the surface, which is more important for higher atomic collision energies. The MS PES favours more than the FPLEPS one the exchange of energy, which can explain its larger adsorption probabilities.
At thermal conditions, the general trends for the probabilities are similar than the ones found for low initial collision energies. FPLEPS and MS probabilities are almost coincident at high temperatures whereas at lower ones the MS adsorption probability decreases when temperature increases due to the physisorption minimum.
The experimental hyperthermal polar scattering angle distribution peak is practically matched for both PES calculations, but the shapes of the theoretical distributions are wider than the experimental one. However, the influence of the experimentally preadsorbed O atoms, which produces O 2 molecules, could affect this final distribution; a direct comparison with the present calculations of atomic scattering over a clean graphite surface could not be totally exact.           The distributions are normalized to unit at the peak.