H 2 Transformations on Graphene Supported Palladium Cluster: DFT-MD Simulations and NEB Calculations

: Molecular dynamics simulations based on density functional theory were employed to investigate the fate of a hydrogen molecule shot with different kinetic energy toward a hydrogenated palladium cluster anchored on the vacant site of a defective graphene sheet. Hits resulting in H 2 adsorption occur until the cluster is fully saturated. The inﬂuence of H content over Pd with respect to atomic hydrogen spillover onto graphene was investigated. Calculated energy barriers of ca. 1.6 eV for H-spillover suggest that the investigated Pd/graphene system is a good candidate for hydrogen storage.


Introduction
Molecular hydrogen is fundamental for a large number of catalytic processes, involving chemical and petrochemical conversions of industrial interest. Included in these, ammonia and methanol syntheses, hydroformylation and Fischer-Tropsch reactions as well as catalytic reforming, hydrotreating (hydrodesulfurization), hydrocracking and either hydrogenation or dehydrogenation of hydrocarbons, oils and fats have to be mentioned [1]. Among the catalytic transformations involving hydrogen, a significant amount of important catalysts contains supported transition metals [2].
A lot of information has been acquired both by experimental and theoretical studies [14][15][16][17][18][19][20][21][22][23][24][25][26], but the analysis of the causes, which dynamically originate the local surface coverage, is still poor [24,27]. In this context, as an example, hydrogen spillover is simply defined as the migration of H atoms to the support from the metal particles, following the dissociative adsorption of H 2 molecule on the metallic sites [28]. However, the spillover mechanism, with the complementary reverse spillover, are still far from being fully understood and further investigations are needed to reach the whole grasp of these phenomena.
From a theoretical point of view, Density Functional Theory (DFT) was employed to compute the energy barrier characterizing the hydrogen migration through the Pt 4 , Pd 6 and Pd 13 clusters' sites and from these to the graphene surface. In all the evaluated cases, the computed energy barriers were 2.0-2.7 eV [29][30][31][32], suggesting that the hydrogen spillover is mostly an energetically hard process. Rangel et al. [33], studying the hydrogen fragmentation, migration and diffusion on one Pd 4 cluster supported on defective graphene, however, observed that the dissociation of the molecular H 2 on the Pd clusters could already occur overcoming energy barriers lower than 0.6 eV and found that the atomic hydrogen migration from the metal to the support occurred overcoming energy barrier values equal to 0.8 eV (on pyridinic defects) and just 0.5 eV (on pyrrolic defects). As a consequence, it was possible to infer that hydrogen fragmentation and migration in these systems might occur almost spontaneously even at room temperature [33].
Studying the energies of sequential H 2 molecules adsorption on Pd clusters at DFT level, Ramos-Castillo et al. [34] suggested that one Pd 4 cluster is able to coordinate several H 2 molecules by shaping H 2 molecular complexes. In this case, the maximum number of H 2 molecules for one Pd 4 cluster supported on graphene resulted equal to four. Conversely, Ebrahimi et al. [35] employed a series of accurate Molecular Dynamics (MD) simulations to characterize the phenomena involved in the H 2 adsorption on graphene nanoribbons, thus observing that the deformation of the graphene structure is closely related to the hydrogen surface coverage.
Recently, Prestianni et al. [24], in an article from now on referred as paper I, introduced a new, quite simple, approach in which the H 2 molecules were addressed "as projectiles" on a supported metal crystallite, namely on a Pd 4 cluster supported on graphene. At lower values of the H 2 kinetic energies, it was observed that the molecular sticking, hence the cluster site hitting was ruled by steering effects imputable to the Pd cluster, which easily reoriented the incoming H 2 molecule to the metallic surface sites, and in this way prompting the fragmentation thus atomic adsorption. Sticking decreased at higher molecular kinetic energies and it appeared that very fast molecules were able to acquire a kind of elusiveness toward the metallic surface which, because of this, was not able to steer and adsorb them.
To add new insights useful to unravel phenomena related to: (i) the surface ruling and hitting of H 2 molecules on metal sites of both pristine and already adsorbed clusters; (ii) the sticking and fragmentation and/or diffusion of atomic hydrogen; and (iii) the existence of spillover and reverse-spillover phenomena, the simple theoretical protocol described in I was mostly employed in the simulations reported in the following.
By these, the individual surface events on a very basic system-taking into consideration the collision of H 2 molecules on a supported fragment, shaped in a periodic system characterized by one unit formed by a (pristine or hydrogen preadsorbed) Pd 4 cluster placed on a C 47 graphene model-were visually evaluated, by employing DFT Molecular Dynamics.
In details, the modeling routine implied the analysis of the local rearrangements following the single collisions of H 2 molecules with a [2 mH+(n−−m)H 2 ]/Pd 4 /C 47 supported fragment, being n + m = 0 − −5. Hence, the effects of the fragment hydrogen coverage as well as the starting speed and the orientation of the H 2 molecule were right correlated to the structural, mechanistic and energetic local changes induced on the same [2 mH+(n−−m)H 2 ]/Pd 4 /C 47 supported fragments.

Results and Discussion
The model in this work is the same already used for the preliminary investigations reported in I. Briefly, the studied system is a tetrahedral Pd 4 cluster anchored on the vancancy of a defective graphene sheet, simulated by a periodic system with a supercell of 47 carbon atoms, with cell vector (12.78, 9.84, 30.00 Å); the projectiles are H 2 molecules having predefined kinetic energies, k e , of 0.005, 0.01, 0.015 and 0.1 eV. In paper I it is evidenced that the hitting of H 2 on Pd 4 /C 47 results in fragmentation only for the lowest values of projectile kinetic energy; the H 2 molecule is adsorbed for k e = 0.015 eV and it bounces at the higher energies. This behavior was ascribed to an energetic balance between the projectile k e and its interaction energy with the cluster. In agreement with the steering mechanism proposed by Gross and Scheffler [36], low kinetic energy projectiles can undergo more easily to a process of steering and reorienting, which could promote their activation and fragmentation. On the other hand, when the second hitting was considered, i.e., H 2 shot either on H 2 Pd 4 /C 47 or on 2HPd 4 /C 47 , it was revealed how the presence of already adsorbed species on the cluster affects the fate of the incoming molecule. In particular, it was evident that, due to repulsive effects from the preadsorbed not-dissociated H 2 , steering effects occur at higher projectile kinetic energies.
The target from which the investigation here reported starts is the system where one H 2 molecule is adsorbed atop the Pd(1) atom (with Pd-H distances of 1.853 and 1.887 Å, and H-H bond length equal to 0.849 Å), one H atom is located on the Pd(1)Pd(3) edge (being 1.887 and 1.738 the H-Pd(1) and H-Pd(3) distances, respectively) and another H atom on the Pd(1)Pd(4) edge (with H-Pd(1) and H-Pd(4) equal to 1.862 and 1.724 Å) of a Pd 4 tetrahedron anchored on defective graphene (see Figure 1 for atoms numbering). To simplify the discussion, the following notation is used throughout: (i) H n 2 is a hydrogen molecule adsorbed on the Pd(n) site of ] site(s) of the cluster; and (iii)the Pd 4 /C 47 cluster/graphene support is implied. With this notation, the target system above is indicated as H 1 2 H 13 H 14 . Clearly, µ 1 , µ 2 and µ 3 hydrogen atoms are characterized by 1, 2 and 3 apices, respectively. Furthermore, it is worth investigating the effect of H 2 adsorption on the Pd 4 configuration with respect to the graphene defect. The C atoms defining the defect contour in C 47 are labeled with letters ranging from a to l (Figure 1), and the anchoring of Pd 4 can be defined by the bonds occurring between the Pd centers of the cluster and those C atoms. Thus, in the H 1 2 H 13 H 14 system, Pd (2) is bonded to the contour atoms e and i (Pd-C = 1.989 and 2.002 Å), Pd(3) is bonded to the contour atom c (Pd-C = 2.383 Å) and Pd(4) to a and l (Pd-C = 2.258 and 2.424 Å); the shorthand notation (2ei3c4al) will be used for this anchoring configuration. For comparison, the one in Pd 4 /C 47 , without adsorbed hydrogen, is (2aei3abc4al), with Pd-C distances equal to 1.937, 2.017, 2.053 Å (2aei), 2.523, 2.396, 2.536 Å (3abc) and 2.027, 2.218 Å (4al), so that in our initial target some bonds between the cluster and the graphene sheet are elongated/broken.  system, of which, in the following, NVT equilibration is analyzed. In this system, the anchoring configuration of Pd 4 on the graphene defect slightly changes, becoming now (2ei3bc4a). At lower energy, 0.005 eV, the H 2 molecule interacts with the cluster but cannot bind to it, hence it steps away. At higher energy, 0.1 eV, the projectile undergoes a sort of reflection: it hits the cluster after 97 fs from the beginning of simulation and its interaction with Pd (1) (3), turns around the cluster and interacts with it in proximity of Pd(4). The palladium cluster, after the hitting, is anchored on graphene only through the Pd(2) and Pd(3) centers. Following equilibration, the H 2 molecule (former the projectile) locates between Pd(1) and Pd(4), and the cluster adopts a (2aei3d4a) position with respect to the defect. Essentially the same result is obtained at k e = 0.01 eV, even if in this case the projectile firstly interacts with Pd (1), and also at k e = 0.015 eV, where following the hit H 2 weakly interacts with Pd(3) and only after equilibration it shifts on the Pd(1)-Pd (4) (1) but one of its atoms interact also with Pd(4). Presumably, this occurrence is observed because both Pd(1) and Pd(4) are already occupied by another H 2 molecule in atop position. Some shapshots of these simulations are collected in Figure 2.   Figure 3). At higher energy (0.01, 0.015 and 0.1 eV), only almost perfect reflections occur. In this case, we also investigated the effect of very high energy particles, namely 0.5 and 1 eV, to see if spillover would take place. Instead, the high energy H 2 always causes the desorption of molecules from the cluster. In particular, the k e = 0.5 eV projectile adsorbs on Pd(1), substituting itself to the two H 2 molecules locating atop it, which at once desorb, while in the k e = 1 eV case the projectile reflects and one H 2 from Pd(1) moves away from the cluster.

Hydrogen Hopping to Graphene
The hopping of one hydrogen atom from the Pd 4 cluster to position C(4) of C 47 (see Figure 1) onto graphene in any direction coherent to our sheet model. The change is legitimated by the fact that the result of hitting by MD simulations on a given H-Pd 4 /C 47 system is to be considered only as one among a multitude of possible outcomes of the real hydrogen hitting, and the adsorption of H 2 on Pd(3) instead of Pd(4) can surely be included amid these outcomes. Further, with this choice, in all of the three systems the H-hopping could occur with somehow the same sequence of steps: the breaking of H-Pd(1) bond, the revolution around Pd(4), the breaking of H-Pd(4) bond and the formation of H-C(4). As a matter of fact, this allows to compare the hopping energetics in the three systems. In the final products the configuration of hydrogen species on the Pd 4 cluster slightly changes, giving rise to the systems H 13  , whose energies show that the spillover reactions always occurs with high endoergicity. Here, the apex C(4) underlines the displacement of the corresponding H atom from the cluster to the graphene site indexed as 4 (see Figure 1), where it binds at a distance of 1.13-1.14 Å.
The results of the NEB calculations from reactants to products are collected in Figure 4. As can be noticed, no matter the starting system, a metastable intermediate locates along the hopping path, where an H atom is located atop the Pd(4) center; the energy of the transition state from the reactant to this intermediate species slightly exceeds 65 kJ mol −1 in the least favorable case, but is always only a few kJ mol −1 above the energy of the intermediate itself. The step from intermediate to product would occur by overcoming an energy barrier of 104.3, 103.2 and 98.4 kJ mol −1 for the paths starting from R1, R2 and R3, respectively. This would suggest that the presence of hydrogen species on the Pd 4 cluster has only a small influence on the H-hopping. On the other hand, the intermediate-to-product step only reflects the energetics of H-Pd(4) bond breaking, while the effect of other H-species on the cluster should be better characterized by taking into account the highest energy barrier with respect to the reactant. If this is done, the difference between the three systems is clearly magnified, indicating that, even if the activation toward H-hopping is always quite difficult, the energy barrier sensibly decreases when the number of H-species on the cluster increases. As a matter of fact, the intermediate can be seen as a small ripple in the minimum energy path from reactant to product and it is possible to think that, in order for the hopping to occur, the reactant has to acquire all in one the energy to overcome the largest barrier.
In the present whole saturated H final product. This last is 140.9 kJ mol −1 less stable than the initial reactant; the higher calculated endoergicity, with respect to the other spillover products, is presumably due to the fact that in this case, if the net process is considered, the H atom hopped on graphene was removed from a stable µ 3 position. Hence, it is possible to argue that for the same reason the energy barrier for this process is the highest of all the investigated ones, which suggests that H-spillover from a fully-saturated palladium cluster could be more difficult due to a necessary rearrangement of the hydrogen atomic species that need to get a right constellation, close to the receiving carbonaceous surface. Finally, NEB calculations were performed to investigate the ease of diffusion of hydrogen atom far away from the anchored cluster, i.e., the hydrogen hopping along the 4→3→2→1 path of graphene (see Figure 1). This was done by considering the spillover product H 1 2 H 134 H C(4) as representative, which has energy of 108.1 kJ mol −1 higher than the R2 reactant. The results, reported in Table 1, show that not all C position in C 47 are equivalent, as could be reasonably expected since the bonding of one H atom on graphene site causes a curvature of the sheet, and the ease with which this occurs depends on the position with respect to the defect. Indeed, in the present case, the properties of the graphene curvature are also heavily affected by the palladium cluster anchored on the defect and by the size of the model used. The same arguments apply to the energies of the transition states for the shifting from one C position to the next, whose trend reflects somehow the stability of the corresponding product and is quantified around 100 ± 15 kJ mol −1 .

Models and Computational Details
Molecular dynamics simulations were performed within the density functional theory framework by using the SIESTA approach as implemented in the homonym code, version 4.2.1 [37]. The employed relativistic norm-conserving Troullier-Martins pseudopotentials and the double zeta plus polarization quality numerical basis sets, on which the SIESTA method is based, were generated and tested as described elsewhere [13]. The PBE exchange-correlation functional in the spin-polarized form was chosen for all calculations, along with a 4 × 4 × 1 sampling of the Monkhorst-Pack grid and a value of 450 Ry for the mesh cutoff.
The algorithm devised to simulate the hitting of H 2 on the target can be summarized as follows. Start with the optimized geometry of Pd 4 /C 47 (the H 1 2 H 13 H 14 system, in the present work) and perform a 10 ps MD simulation (timestep of 0.5 fs) in the NVT (298 K) ensemble. (A) Once the system is equilibrated, take the atomic positions, velocities and forces of the last snapshot and add to the system a H 2 molecule with the following characteristics: the first H atom is placed at a distance of 7 Å from the center of mass (CM) of the palladium cluster, with a given angle with respect to the normal to surface; the second H atom is located in the same direction at a distance of 0.75 Å from the first; the velocity vectors of each H atom have magnitudes according to the initial kinetic energy of the projectile and are oriented toward the CM of Pd 4 ; and, finally, the forces acting on the two H centers can be safely assumed to be negligible (hence set to zero) due to the distance of H 2 from the target. Restart the MD simulation for a duration of 2 ps, now using the NVE ensemble in order to not redistribute the energy possessed by the projectile, and analyze the product of the hitting. Do this for projectiles having initial kinetic energies of 0.005, 0.01, 0.015 and 0.1 eV. Since the aim is to investigate the adsorption and/or the fragmentation of H 2 on Pd 4 /C 47 , the simulations ended with simple reflection or bouncing of the projectile were discarded, while those ended with an effective hitting, i.e., an incipient sticking, were retained, analyzed and then subjected to equilibration in the NVT ensemble for 10 ps. The same procedure is repeated starting from (A) to simulate the hitting of a new H 2 molecule and continue to shoot H 2 until no effective hit occurs, i.e., the Pd 4 cluster appears to be H-saturated.
Hopping of hydrogen species from the metal cluster to the graphene surface sites (spillover) was studied by means of a home-made code, which implements the Climbing Image Nudged Elastic Band method (CI-NEB) [38,39] and is interfaced with the SIESTA program. Reactants and products were fully optimized by SIESTA periodic calculations and NEB images were generated using the image depending pair potential [40]. Every NEB calculation used 13 images connected each other by dynamic springs and a convergence threshold of 3 × 10 −3 E h /Å on the norm of atomic forces was employed.

Conclusions
DFT-based molecular dynamics simulations show that a small palladium cluster anchored on a defective site of graphene can be saturated with hydrogen by means of a flux of low-energy H 2 molecules. Cluster H-saturation is reached when five hydrogen molecules are adsorbed/fragmented, taking into account that some cluster sites are not available for the interaction with H 2 . In fact, even if the cluster anchoring configuration slightly changes and/or weakens according to the degree of H-adsorption, one face of the cluster is always firmly interacting with the graphene defect region. As the number of hydrogen species on palladium increases, the energy barrier that rules the H-spillover onto graphene slightly decreases. H-spillover remains an onerous process, which is not surprising given the high endoergicity of the reaction, but the present results indicate that, in the here investigated system, a small palladium cluster on defective graphene, the energy barriers are reduced by 0.5-1.0 eV with respect to other systems (excluding of course those where H-spillover occurs toward a graphene defect). This appears true also if the spillover requires some sort of redistribution/rearrangement of atomic hydrogen on the cluster.
Finally, by considering the whole web of the different pathways presented in this study, straightforwardly it emerges that the reverse-spillover (i.e., hopping from graphene to Pd 4 , rearrangement of hydrogen species on the cluster, and formation and desorption of H 2 ) is supposed to be an easy process, when occurring on the graphene-anchored palladium clusters. This should encourage experimental investigations to deeper probe this kind of systems, hence to give them a chance as reversible hydrogen storage device.

Conflicts of Interest:
The authors declare no conflict of interest.