The mechanism of chemisorption of hydrogen atom on graphene: Insights from the reaction force and reaction electronic flux Diego Cortés-Arriagada, Soledad Gutiérrez-Oliva, Bárbara Herrera, Karla Soto, and Alejandro Toro-Labbé Citation: The Journal of Chemical Physics 141, 134701 (2014); doi: 10.1063/1.4896611 View online: http://dx.doi.org/10.1063/1.4896611 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/13?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Adsorption of nitrogen oxides on graphene and graphene oxides: Insights from density functional calculations J. Chem. Phys. 134, 044710 (2011); 10.1063/1.3541249 Quantum study of Eley-Rideal reaction and collision induced desorption of hydrogen atoms on a graphite surface. I. H-chemisorbed case J. Chem. Phys. 124, 124702 (2006); 10.1063/1.2177654 The importance of hydrogen’s potential-energy surface and the strength of the forming R – H bond in surface hydrogenation reactions J. Chem. Phys. 124, 044705 (2006); 10.1063/1.2159482 Dissociative adsorption of NO upon Al(111): Orientation dependent charge transfer and chemisorption reaction dynamics J. Chem. Phys. 117, 8185 (2002); 10.1063/1.1519107 Atomic and molecular hydrogen interacting with Pt(111) J. Chem. Phys. 111, 11155 (1999); 10.1063/1.480473

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00

THE JOURNAL OF CHEMICAL PHYSICS 141, 134701 (2014)

The mechanism of chemisorption of hydrogen atom on graphene: Insights from the reaction force and reaction electronic flux Diego Cortés-Arriagada,a) Soledad Gutiérrez-Oliva, Bárbara Herrera, Karla Soto, and Alejandro Toro-Labbé Nucleus Millennium Chemical Processes and Catalysis, Laboratorio de Química Teórica Computacional (QTC), Departamento de Química-Física, Facultad de Química, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, Macul, Santiago, Chile

(Received 19 June 2014; accepted 16 September 2014; published online 1 October 2014) At the PBE-D3/cc-pVDZ level of theory, the hydrogen chemisorption on graphene was analyzed using the reaction force and reaction electronic flux (REF) theories in combination with electron population analysis. It was found that chemisorption energy barrier is mainly dominated by structural work (∼73%) associated to the substrate reconstruction whereas the electronic work is the greatest contribution of the reverse energy barrier (∼67%) in the desorption process. Moreover, REF shows that hydrogen chemisorption is driven by charge transfer processes through four electronic events taking place as H approaches the adsorbent surface: (a) intramolecular charge transfer in the adsorbent surface; (b) surface reconstruction; (c) substrate magnetization and adsorbent carbon atom develops a sp3 hybridization to form the σ C-H bond; and (d) spontaneous intermolecular charge transfer to reach the final chemisorbed state. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4896611] I. INTRODUCTION

The adsorption of atoms and molecular system on surfaces is today a broad field of study mainly in the effort to understand the processes whereby adsorbates interacts on the adsorbents, chemical modification of surfaces allowing adsorption of new analites, or mechanism whereby adsorbates form new molecules. In this sense, a considerable amount of experimental and theoretical studies has been developed to analyze the interaction of hydrogen (and hydrogen molecule) onto graphitic surfaces. The aim is mainly from its implications in use for hydrogen storage,1, 2 study of hydrogen-wall interactions in reactors for nuclear fusion technology,3–5 band gap tuning for nanoscale electronic devices,6, 7 and to clarify the formation of hydrogen molecule in the interstellar media proposed to take place by recombination on graphitic dust grains.8–11 As noted above, several theoretical studies have been developed to characterize the interaction of hydrogen onto graphitic surfaces employing polycyclic aromatic hydrocarbons (PAHs) (as coronene or pyrene) and periodic surfaces of graphite or graphene.5, 8, 10, 12–19 Mainly, a lot of density functional theory (DFT) methods were employed in order to account for accurate values of energetic parameters, which has shown to be in good agreement with rigorous coupled cluster methods like CCSD(T)/cc-pVTZ.15 From these studies, it was predicted that H physisorption onto graphitic surfaces takes place at ∼3.0 Å onto the surface. For chemisorption, the hydrogen atom need to overcome a predicted barrier of ∼4.6 Kcal/mol, confirmed from experiments by high-resolution electron energy loss spectroscopy,20 forming a transition state onto a C atom (so-called top site) at a) Author to whom correspondence should be addressed. Electronic mail:

[email protected] 0021-9606/2014/141(13)/134701/9/$30.00

∼1.7–1.8 Å. In the chemisorbed state, hydrogen atom bonds to C atom and in the process the latter develops sp3 hybridization and puckers out of the surface basal plane; therefore, the substrate reconstruction is necessary to allow the hydrogen atom chemisorption on graphene.8, 21 All process has been predicted with a reaction energy of ∼12.5 Kcal/mol15 and generates on graphene the midgap state due to a pz unpaired electron.22 All of these analyses have focused into study only the stationary points on the potential energy surface, being necessary an extended analysis along the reaction coordinate to gain insights about the overall processes taking place during the adsorption. Expecting the adsorption on graphene involves processes such as charge transfer, bond weakening, and bond strengthening, we can characterize them from chemical descriptors as those proposed from the DFT. In recent years, reaction force23 and reaction electronic flux24, 25 has been used to rationalise mechanisms involved in a wide type of chemical reactions as SN 2, proton transfers, isomerisations, among others.26–28 However, its use in adsorption problems involving physical and chemical interactions has not been explored yet. In this work, by employing a dispersion corrected DFT method, we use the reaction force and reaction electronic flux to study the mechanism whereby hydrogen is chemisorbed onto graphitic surfaces. In addition, to account for the reliability of the chemical information obtained by these descriptors, we compared it to that coming from electron population analysis. II. THEORY A. Reaction force

In the course of a chemical or physical process, the force F is defined as the negative of the energy derivative along the

141, 134701-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00

134701-2

Cortés-Arriagada et al.

J. Chem. Phys. 141, 134701 (2014)

until reaching the final state. In summary, although structural and electronic effects are present all along the process R → P, it has been shown that structural effects are dominant in the reactant and products regions whereas changes in electronic properties dominate within the transition state region. Integration of F(ξ ) within the reaction regions allows obtaining a decomposition of the reaction (E◦ ) and activation (E= ) energies. In the case of a single step reaction, one can define: ξ1 W1 = −

ξT S F (ξ )dξ > 0,

ξ2

ξP F (ξ )dξ < 0,

ξT S

reaction coordinate ξ :23 dE . F (ξ ) = − dξ

(1)

F(ξ ) is refereed as the reaction force; this global property is zero at every maximum or minimum of the potential energy profile (E(ξ )) and presents critical points where the derivative of the E(ξ ) is zero. In a single step reaction R → P, where reactants (R) converts into products (P) passing through a transition state (TS), two critical points are found: at ξ 1 between ξ R and ξ TS and at ξ 2 between ξ TS and ξ P (Fig. 1). By means of ξ 1 and ξ 2 it is possible to define the reaction regions in which different phenomena might be taking place; these are the reactant [ξ R , ξ 1 ], transition state [ξ 1 , ξ 2 ], and product [ξ 2 , ξ R ] regions. The physical interpretation of F(ξ ) emerges through the analysis of the mechanism that might be driving the reaction at every reaction region. At the reactant region, in the [ξ R , ξ 1 ] range there is a retarding force reaching a maximum value at ξ 1 . In the transition state region [ξ 1 , ξ 2 ], appears a driving force balancing the retarding force and reaching a maximum at ξ 2 , passing by the TS where retarding and driving forces are cancelled at ξ TS . In the final step, between [ξ 2 , ξ R ], the driving force decreases as the products are obtained. Studies of many kinds of chemical reactions26–32 have shown the main effects dominating in each region. Reactant region is preparative in nature and structural effects associated with geometrical changes in the chemical species (i.e., bond stretching, angle bending, geometrical changes due to atomic re-hybridization, etc.) prevail, in agreement with the retarding character of the F(ξ ) in this region. In the transition state region, the driving force is mainly associated with electronic reordering (coming from bond formation/strengthening, bond break/weakening, charge transfer, polarization, etc.). Most of the electronic activity that occurs in a chemical reaction takes place in the transition state region. Finally, in the product region, structural relaxation takes over and drives the reaction

F (ξ )dξ > 0, (2) ξ1

W3 = − FIG. 1. Potential energy (solid line) and reaction force (dashed line) profiles of an exergonic single step reaction as a function of the reaction coordinate ξ . On top are displayed the boundaries for the reactant (R), transition state (TS), and product (P) regions as defined by the minima and critical points in the energy and reaction force profiles, respectively.

W2 = −

ξR

W4 = −

F (ξ )dξ < 0, (3) ξ2

and therefore E◦ and E= of a single step reaction can be expressed as   E ◦ = E(ξP ) − E(ξR ) = W1 + W2 + W3 + W4 , (4)   E = = E(ξT S ) − E(ξR ) = W1 + W2 .

(5)

From Eqs. (4) and (5) it is possible to analyze the predominant effects on E◦ and E= from the relative weight of the involved reaction works and character of each region. In this context, W1 and W4 are mainly structural energies while W2 and W3 are mainly associated to electronic processes. For instance, the energies required for the structural changes in reactant region (W1 ) plus energy mainly associated to electronic reordering in the first stage in the TS region (W2 ) define the activation energy, and the activation energy appears to be defined through a balance between structural and electronic effects. In this context the nature of the activation energy can be characterized through the relative contribution of W1 and W2 to E= . B. Reaction electronic flux

Conceptual DFT provides a set of properties that characterize both global and local reactivity of molecules, they are obtained from the derivatives of the energy E with respect to the number of electrons N and the external potential υ(r).33 At the global level, the electronic chemical potential characterizes the response of E with respect to changes in N, at constant υ(r). Chemical potential (μ) is related to the electronegativity χ thought μ = −χ .34   ∂E ∼ μ= (6) = −χ . ∂N υ(r) Use of the finite difference approximation allows to obtain the chemical potential in terms of the first ionization potential (IP) and the electron affinity (EA). Further approximation using the Koopmans Theorem leads to μ in terms of the

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00

134701-3

Cortés-Arriagada et al.

J. Chem. Phys. 141, 134701 (2014)

eigenvalues ε of the frontier orbitals [HOMO (H) and LUMO (L)]:   1 ∂E 1 μ= ≈ − (I P + EA) ≈ (εH + εL ). (7) ∂N υ(r) 2 2 Note that Eq. (7) is only applicable for closed-shell systems. For the open-shell case we need to use a general form as defined from the spin polarized DFT:35   β β α + εH 1 εH εLα + εL ∼ μ= + , (8) 2 2 2 where α and β stand for the spin state. In the closed-shell case β β α = εH and εLα = εL , therefore Eq. (8) is reduced to Eq. (7). εH In this context, profiles of chemical potential along a reaction coordinate might be obtained quite easier. The reaction electronic flux determines the changes of μ along the reaction coordinate ξ :24, 25 J (ξ ) = −

dμ . dξ

(9)

It has been shown that J(ξ ) is related with the spontaneity of the electronic activity taken place along the reaction coordinate. In this sense, positive values of REF indicate spontaneous changes in the electronic density driven by bond strengthening or bond formation processes; while, negative values of REF indicate that non-spontaneous changes of the electron density driven by bond weakening or breaking processes. In addition, although charge transfer and polarization are phenomena difficult to separate during a physical or chemical process, a partition has been proposed25 in the context of REF theory as follows: J (ξ ) = Jp (ξ ) + Jt (ξ ),

(10)

where [Jp (ξ )] and [Jt (ξ )] stand for polarization and transfer contributions. A partition of the molecular system in i molecular fragments allows obtaining Jp (ξ ) and Jt (ξ ) by fragment contributions:25 Jp (ξ ) =

i=n i=1

Jpi (ξ )

 i=n  Ni dμi , =− N dξ i=1

(11)

Jt (ξ ) =

i=n i=1

Jti (ξ )

 i=n   Ni d  μi (ξ ) − μ(ξ ) , = N dξ i=1

(12)

where Ni and N are the number of electrons at fragment i and in the whole system, respectively; μi stands for the chemical potential of fragment i. The condensed polarization flux Jpi (ξ ) is related to the electron density deformation of fragment i caused by external field of the other fragments; Jti (ξ ) accounts for inter-fragments electron transfer processes. In order to obtain μi , frontier orbital eigenvalues of the fragment i with the basis set of the whole system are used, as obtained by the counterpoise method.36 An extended revision of REF and reaction force theories can be found in the original works.23–25 III. COMPUTATIONAL DETAILS

For all simulations, and in order to account for a model that allows an adequate description of reaction force and reaction electronic flux profiles, we adopted three types of polycyclic aromatic hydrocarbons (PAHs) as finite models for graphene (Fig. 2): pyrene (C16 H10 ), coronene (C24 H12 ), and a zigzag model called G4×4 (C48 H18 ). At the DFT level, the PBE37 functional was used in combination with the Dunning correlation consistent polarized double zeta basis set (cc-pVDZ). PBE functional was selected considering its performance on non-covalent interactions of adsorbates upon graphene in reported and our previous works.38, 39 The spin unrestricted formalism (UPBE) was used for calculations involving open-shell electronic states. In these cases, mixture of higher spin states was found to be negligible; the maximum obtained value for the spin-squared operator S(S + 1)(where S is the spin angular momentum) was 0.76 in comparison to an expected value of 0.75, discarding spin contamination problems during calculations (Fig. S1 in the supplementary material).52 Dispersion force effects for energy and gradients were included by the pairwise dispersion correction DFT-D3, in combination with the Becke-Johnson damping scheme in order to avoid repulsive interatomic forces at short distances.40–42 Relaxed potential energy surface (PES) scans were performed along of the reaction coordinate C1 H1 , obtaining 60 points from 1.0 to 3.5 Å. In geometry optimizations, the Broyden–Fletcher–Goldfarb–Shanno Hessian update method43 was adopted. In geometry optimization and

FIG. 2. PAHs used on all calculations. Atomic hydrogen is adsorbed onto C1 . Hydrogen atoms at the edges were deleted in order to simplify.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00

134701-4

Cortés-Arriagada et al.

J. Chem. Phys. 141, 134701 (2014)

FIG. 3. Hydrogen interaction on all graphene models as an adsorption process, depicting reactants (physisorbed state), transition state (TS), and products (chemisorbed state). Distances are in angstroms (Å).

SCF steps, convergence tolerance values of 1 × 10−6 and 1 × 10−8 hartree were used, respectively. Vibrational frequencies were obtained for stationary points in order to ensure for the nature of reactants, products, and transition states. Atomic net charges, Wiberg bond orders, and hybridization indexes were obtained from the Natural Bond Orbital (NBO) Analysis. Basis set superposition errors (BSSE) were corrected by the counterpoise (CP) method by Boys and Bernardi.36 Structure optimizations were carried-out in the ORCA 3.0.1 program,44 and properties calculations were performed in the Gaussian 09 program.45

IV. RESULTS AND DISCUSSION A. Potential energy and reaction force profiles

The results consider the adsorption processes as a function of the vertical distance C1 -H1 . Fig. 3 shows the hydrogen

interaction with all graphene models as an adsorption process, depicting reactants (physisorbed state), transition state (TS), and products (chemisorbed state). Fig. 4 shows the potential energy and reaction force profiles and Table I shows the reaction and energetic parameters for all analyzed models. Note that the process under study is physisorption → chemisorption and the reaction goes from higher to lower values of the C1 -H distance. In this context, the works defined through Eqs. (2) and (3) change signs. Firstly, H physisorption take place at ∼3.0 Å on C1 with an exergonic energy (Ephys ) of −0.68 to −0.80 Kcal/mol, in agreement with a predicted value of −0.7 Kcal/mol from DFT calculations using the BP86-D/def2-TZVP method17 ; note that the value from the G4×4 model shows the best agreement with experimental value of −0.90 Kcal/mol for H physisorption on graphite.46 Although both top and hollow sites have been theoretically determined to be the minimum in the physisorbed state, disposition of H is predicted

FIG. 4. (a) Energy and (b) reaction force profiles in Kcal/mol for hydrogen interaction on all PAHs models. Vertical lines in reaction force account for the limits among reactants, transitions state, and product regions. W represents the reaction works.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00

134701-5

Cortés-Arriagada et al.

J. Chem. Phys. 141, 134701 (2014) =

TABLE I. Physisorption energy (Ephys ), reaction energy (E◦ ), activation energy (E= ), reverse activation energy (Er ), and reaction works (Wn ) according to reaction force theory. All values include counterpoise and dispersion force correction and are expressed in Kcal/mol. Ephys = E(PAH − H)phys − (EPAH + EH ); E ◦ = E(P AH −H )chemi − (EP AH + EH ) = Ephys − (W1 + W2 + W3 + W4 ); E = = E(P AH −H )ts =

− E(P AH −H )phys = −(W1 + W2 ); Er = E(P AH −H )ts − E(P AH −H )chemi = W3 + W4 ; where E(PAH-H)phys , E(PAH-H)chemi , and E(PAH-H)ts are the energies of the PAH-H system in the physisorbed, chemisorbed and transition state, respectively. EPAH and EH are the energy of the isolated PAH and hydrogen atom, respectively. Adsorbent

Ephys

E◦

E=

E= r

W1

W2

W3

W4

Pyrene Coronene G4×4 Experimental

−0.68 −0.73 −0.80 −0.9a

−13.88 −14.89 −19.15 −10.3b,c

5.25 5.17 4.64 4.6d

18.45 19.33 22.99 14.0c,e

−3.78 (72%) −3.77 (73%) −3.34 (72%) −3.3f

−1.47 (28%) −1.40 (27%) −1.30 (28%) −1.3f

11.68 (63%) 12.38 (64%) 15.40 (67%) 9.4f

6.77 (37%) 6.95 (36%) 7.58 (33%) 4.6f

a

Experimental values were obtained from Ref. 46. Estimated from Refs. 20, 46, and 47. Note that inclusion of ZPE corrections affects energy parameters in ∼4 Kcal/mol when chemisorbed state is involved. d Experimental values were obtained from Ref. 20. e Experimental values were obtained from Ref. 47. f Estimated from percentage associated to reaction works of the H-G4×4 system. b c

to be in the top site on C1 atom, which is even a discussion focus from DFT and ab initio calculations.10, 15–17 In the next stage, H chemisorption on top site of all adsorbents leads to exergonic products over-passing a chemisorption barrier that ranges from 5.25 to 4.64 Kcal/mol, depending on the adsorbent and in good agreement with the experimental value of ∼4.6 Kcal/mol.20 The smallest barrier is obtained with the G4×4 model, indicating the importance of dispersion forces in an extended surface to allow the chemisorption on graphitic surfaces from the physisorbed state. Besides, this barrier is high enough to prevent chemisorption at ambient conditions.14 Transition state is located at a C1 H1 distance of ∼1.8 Å, which agree with distances determined theoretically (1.72–1.80 Å)10, 12, 15–17 ; at this point, C1 atom is placed over the plane of the PAH due to the hybridization change from sp2 to sp3 to form a σ bond between 1s (H1 ) and 2pz (C1 ) orbitals in the chemisorbed state. Finally, chemisorption is reached at ∼1.13 Å in the reaction coordinate, with exergonic chemisorption energies of −13.88, −14.89 and −19.15 kcal/mol for pyrene, coronene and G4×4 , respectively, in agreement with values reported for related models.5, 8, 10, 12–19, 48 A deeper minimum is found for the G4×4 model in ∼5 Kcal/mol compared to pyrene and coronene, indicating the importance of lattice size to determine accurate values of E◦ to compare with graphite or graphene. Indeed, G4×4 chemisorption energy is very close to values determined from extended models as C42 H16 (−19.14 Kcal/mol, B3LYP/6-31G(d,p)12 and graphene under periodic boundary conditions avoiding BSSE (∼−19 Kcal/mol; PBE/plane waves).14, 48 Moreover, although three different sites are commonly proposed for adsorption of atoms on graphene (top, adsorption on carbon atom; bridge, adsorption on C-C bonding; and hollow, adsorption at the centre of benzene type ring),12, 16 theoretical studies show that the adsorption of atoms with one valence electron take place on the top site,49 therefore our results agree with these findings. The reverse energy barrier characterizes the desorption process, which needs over-passing a barrier larger than 18 Kcal/mol, naturally due to that desorption is a bond breaking = process; corrected Er by zero point energy (ZPE) have reasonable agreement with the experimental desorption

activation energy of H bonded on graphite (∼14 Kcal/mol).47 It is important to note that inclusion of ZPE corrections (Table S1 in the supplementary material)52 affects energy parameters in ∼4 Kcal/mol when chemisorbed state is involved, in agreement with thermochemical calculations, which is due to strong H bonding to the surface in the chemisorbed state.17 Reaction force allows fragmenting the reaction coordinate in three regions (reactants, transition state and products) according to critical points located at ξ ≈ 2.0 and ξ ≈ 1.3 Å in the reaction force profile (Fig. 4(b)). By using Eqs. (2) and (3), the reaction works involved at each region along the reaction coordinate are obtained. From W1 and W2 , the works involved in the chemisorption energy barrier (E= ) were found, which would be mainly dominated by structural rearrangements in a 72%–73%. The magnitude of W1 is well related to the substrate reconstruction needed to allow the chemisorption occurring above the TS distance (ξ > 1.8 Å),8 indicating that reconstruction of the host material is the main effect involved in the adsorption barrier, which has been shown to be needed to allow the chemisorption of hydrogen.8, 21 Assuming the same relative weight of W1 and W2 to define E= in the H-G4×4 system, the experimental value of reaction works can be estimated, they are also displayed in Table I. Until TS, adsorption is mainly dominated by dispersion forces, therefore the work associated to electronic effects is lesser than the structural one, with contributions of 27%–28%. After TS, chemisorption is mainly driven by spontaneous electronic rearrangements, with W3 of 64%–67% of the reverse energy = barrier (Er ). Finally, W4 shows structural rearrangements = corresponding to the 33%–36% of Er , respectively; 3% lesser of G4×4 than coronene and pyrene would be due to the fact that the later ones are distorted at the edges near to C1 . Note that W4 and W3 are negatives for the chemisorption process; therefore positive values of W4 and W3 represent the most important works involved in the endergonic H desorption. In summary, reaction force indicates that for the chemisorption process, the structural work is the most important contribution to the energy barrier (W1 > W2 ), while that in the desorption process, electronic work is the greatest contribution (W3 > W4 ).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00

134701-6

Cortés-Arriagada et al.

FIG. 5. Evolution along the reaction coordinate for C1 -CN (N = 2, 3, and 4) distance and C2 -C1 -C3 -C5 dihedral angle. Atom labels are depicted in Fig. 2. Results for G4×4 are presented as representative for all adsorbents (all results included in the supplementary material).52

In order to account for the structural nature of W1 and W4 , we analyze the structural changes along the reaction coordinate associated to the substrate reconstruction (Fig. 5). In addition to C1 -H1 distance, we account for C1 -CN distance, where N is the label for adjacent atoms to C1 (Fig. 2). It can be seen C1 -CN distance starts at 1.42–1.43 Å, according to C-C bond in graphene. From reactants to the first region of TS, C1 -CN distance changes ∼0.02 Å, according with the low barrier for chemisorption. Furthermore, when TS is reached, C1 -CN distance increase dramatically until chemisorbed state is reached; at this point on the surface is found a C1 -CN bond of ∼1.52 Å, near to the value of a diamond type structure. Besides, the dihedral angle C2 -C1 -C3 -C5 measures as C1 atom puckers out of the basal plane changing the surface planarity; this structural change is involved in E= because it starts to decrease at ξ ≈ 2.5Å from 180◦ to ∼160◦ in the chemisorption region where C1 develops sp3 hybridization. The results indicate that both W2 and W3 involve structural components, but electronic effects in TS region will be dominant as discussed below by means of the REF analysis. The above discussion points out that most of the structural effects in W1 are due to reconstruction at the graphene surface. According to the surface reconstruction that dominates the energy barrier of H chemisorption on the graphene models, it could be expected that this barrier will be affected by the doping level. The activation energy is expected to be reduced through controlling the substrate reconstruction before adsorption. This can be done by doping the surface with dopants atoms with larger atomic radii compared to a carbon atom, which avoid surface reconstruction due to doping process is itself highly structural and changes the local structure of the surface. For instance, the energy profiles of H adsorbed on A-doped G4×4 model (where A is Al, Si, or P) show that the adsorption proceeds without an activation energy barrier associated. In this regard, Fukushima et al. report a decrease of ∼3 Kcal/mol in the activation energy of H adsorption on Al-doped graphene50 with respect to the undoped surface. On the other hand, doping with atoms of similar radii compared to carbon, as nitrogen or boron, does not affect the topological structure of the graphitic surface, so that the

J. Chem. Phys. 141, 134701 (2014)

surface must develops reconstruction to allow the chemisorption resulting in an activation energy barrier. In fact, in the case of N-doped G4×4 , we estimate activation energy barrier > 10 Kcal/mol. Furthermore, the decrease of the activation barrier would result in an increase of the binding energy.32 In other respects, Huang et al.51 establish that an increase in the desorption energy barrier is reached by charge doping, which decrease the energy cost necessary to break the aromatic π bonds in graphene and enhancing the strength of the C-H bond in the chemisorbed state. From the chemical point view, the transition state would be more difficult to reach. As shown above, desorption energy barrier is dominated by electronic work, hence for H desorption from a charge doped substrate, it is expected that electronic effects associated with bond break/weakening are increased as the transition state region is increased. Taking into account hole and electron doping in the H-G4×4 system, it was found that W3 increase from 67% (in the undoped substrate) to 80% (in both charge-doped substrates).

B. Reaction electronic flux

Fig. 6(a) depicts the J(ξ ) profiles for H adsorption on all the adsorbents and shows four critical points with the electronic events occurring during the chemisorption. In the reactant region we observe a slightly J(ξ ) < 0 indicating a non-spontaneous electronic activity associated with the bond weakening of C1 -CN bonds; intensity of REF in this region is low due to structural rearrangements and dispersion forces dominate the chemisorption barrier as indicated by the reaction force. Entering to the TS region, two maxima (J(ξ ) > 0) are found: the first one in the TS distance indicating the C1 H1 bond strengthening, and the second one entering into the product region that can be associated to the final electronic arrangement on all the system to reach the final chemisorbed state. It is necessary to point out that from REF theory it can be expected that most of the electronic activity be associated to charge transfer contributions. The polarization effect of H atom on the adsorbent electron density is minimal due to its reduced electron density compared to the adsorbent. Indeed, polarization flux has a minimal activity along the reaction coordinate as shown in Fig. 6(b); therefore, the chemisorption process is dominated by charge transfer effects. On the other hand, according to Eqs. (11) and (12), as N  NH on these systems (ZH = 1), both polarization and transfer effects are dominated by the adsorbent (decomposition by fragment is not necessary). Taking this into account, we can characterize the four electronic events mentioned above by analyzing other independent electronic properties such as atomic charges, bond orders and hybridization degree of the adsorbent carbon atom (Fig. 7). These properties were selected for the G4×4 model as representative for all the adsorbents (all systems are included as the supplemental material).52 From the zero flux regime, chemisorption begins with non-spontaneous [Jt (ξ ) < 0] intramolecular charge transfer from the C1 atom to the adsorbent thus weakening the C1 -CN bonds. This fact is in agreement to increase in the total dipole moment near

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00

134701-7

Cortés-Arriagada et al.

J. Chem. Phys. 141, 134701 (2014)

FIG. 6. (a) Reaction Electronic Flux (REF, J(ξ )) profiles in Kcal/mol Å, for hydrogen adsorption on all adsorbents. (b) Decomposed J(ξ ) profile in transfer (Jt ) and polarization (Jp ) contributions for the H-G4×4 model.

to the TS region, which is due to charge decrease of C1 atom and not H1 or the surrounding C atoms, as noted from charge profile in Fig. 7. This step in the reaction is allowing the substrate reconstruction and sp3 hybridization of C1 in a later step. The two peaks in the Jt (ξ ) into the TS region indicate two spontaneous charge transfer events [Jt (ξ ) > 0]. The first one at the TS distance (∼1.80 Å), mainly associated with the intramolecular electronic reordering due to the surface reconstruction and formation of the TS because of the charge transfer from hydrogen is 0 is in the limit of TS-P regions (∼1.34 Å) coming from the final electronic reordering before to reach the chemisorbed state including charge transfer from the adsorbed H atom (∼0.3 electrons). From the dipole moment, we found local minima at 1.34 Å (0.1 D) and a peak at ∼1.80 Å (0.6 D) indicating important charge density distributions at these points in agreement with the peaks in Jt (ξ ), which are qualitatively related with intensity in the changes of dipole moment. Finally, a poorly defined maximum of pyrene is noted in the TS region of J(ξ ) profile, most possibly due to the system size since coronene and G4×4 exhibit a quite similar REF profile; so we recommend using at least coronene in order to

FIG. 7. Charge of H1 and CN (N = 1, 2, 3, and 4) atoms and total dipole moment in debye (D) along the reaction coordinate for the H-G4×4 system. All systems are included in the supplementary material.52

produce reliable REF profiles of the system. In addition, REF indicates that most of the electronic events are occurring in the TS region, as indicated by the reaction force, and supporting W2 and W3 are mainly electronic works even if structural changes are also taking place in the TS region. C. Bond orders and C1 atom hybridization

Fig. 8 shows the Wiberg bond order derivatives and hybridization of C1 atom along the reaction coordinate for the G4×4 system. The highest intensity of bond order reordering takes place in the TS region (∼1.6 Å); note that weakening (negative derivative) and strengthening (positive derivative) of C1 -CN (parameter associated to the surface reconstruction) and C1 -H1 bonds, respectively, are occurring at the same time, i.e., the process is synchronous. Moreover, from the C1 hybridization profile, we can see the snappish change of sp2 to sp3 occurs at ∼1.6 Å, after the TS, the same point where the bond order derivative reaches its maximum value. These changes along the reaction coordinate match with the minima of J(ξ ) and Jt (ξ ) profiles in the TS region at ∼1.6 Å. This result indicates that at this point occurs the forming of C1 H1 bond, and C1 develops a sp3 hybridization, just after the

FIG. 8. Wiberg bond order derivative of C1 -H1 and C1 -CN (N = 2, 3, and 4) bonds and hybridization of C1 atom along the reaction coordinate for the H-G4×4 system. All systems are included in the supplementary material.52

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00

134701-8

Cortés-Arriagada et al.

J. Chem. Phys. 141, 134701 (2014)

V. CONCLUSIONS

FIG. 9. Natural bond orbitals (NBOs) describing the σ bond at the chemisorbed state. No NBO is found for the C1 -H1 bond at the TS distance, only natural atomic orbitals (NAOs) corresponding to 2pz of C1 atom and 1s of H1 atom.

transition state. Indeed, Fig. 9 shows that there is no NBO associated with the C1 -H1 bond in the TS. After this point, adsorption (or desorption) can go forward by spontaneous charge transfer processes as discussed above. In addition, this information from the REF profile indicates that in the chemisorption process, graphene loses aromaticity after the TS is established, indicating that TS is similar to the reactants. Finally, from Fig. 10, it is interesting to note that the high intensity in the spin density change takes place at the same point discussed above, suggesting the midgap state of graphene or the substrate magnetization, induced by pz vacancy unpaired electron on the adsorbent surface,22 is created simultaneously with the bond weakening/strengthening and where the C1 hybridization change reaches the maximum intensity, as indicated by REF at ξ ≈ 1.6 Å. These last findings allow establish the main process involved in the reverse activation energy, which was found to be dominated by electronic effects. Note that in the desorption process, the events taking place are opposite to those discussed, i.e., change in the spin state of the adsorbent from high to low spin, charge transfer to reach the transition state, and change in hybridization of C1 from sp3 to sp2 .

At the PB3-D3/cc-pVDZ level of theory, the steps involved in the chemisorption of hydrogen onto graphene were fully analyzed from the reaction force and reaction electronic flux. As models of graphene we used pyrene, coronene and a G4×4 cell. Our results indicate that the minimum size model to study the adsorption of hydrogen onto graphene can be represented by coronene. However, an extended surface as G4×4 model is necessary to determine accurate values of adsorption energies and comparing with graphitic surfaces. Reaction energies and barriers were in good agreement with experimental data. Through the reaction force analysis, it was determined that chemisorption energy barrier is mainly dominated by structural work (73%) associated to the substrate reconstruction, while the electronic work is the greatest contribution (67%) to desorption activation energy. On the other hand, reaction electronic flux in combination with the NBO electron population analysis shows that four steps are involved in all chemisorption process: (a) intramolecular charge redistribution in the adsorbent; (b) surface reconstruction and spontaneous charge density distribution; (c) change in hybridization of C1 to allow formation of the C1 -H1 bond, and substrate magnetization (midgap state of graphene) induced by unpaired electron on the adsorbent surface; and (d) final spontaneous intermolecular charge redistribution to reach the final chemisorbed state. Note that electronic activity taking place in steps (c) and (d) explains why the desorption energy barrier is mainly dominated by electronic effects. In this way, reaction force and REF show to be useful tools to study the processes whereby adsorbates interact onto surfaces.

ACKNOWLEDGMENTS

The authors thank FONDECYT through Project Nos. 1130072, 1120093, and 1141098 and ICM Grant No. 120082 for financial support. D.C.-A. acknowledges Postdoctoral FONDECYT Project No. 3140314 and to QTC members R. Duran, S. Miranda, D. Ortega, and N. Villegas for useful discussions. 1 V.

FIG. 10. Spin density derivative corresponding to the chemisorbed H atom (H1 ).

Tozzini and V. Pellegrini, Phys. Chem. Chem. Phys. 15(1), 80–89 (2013). 2 K. Spyrou, D. Gournis, and P. Rudolf, ECS J. Solid State Sci. Technol. 2(10), M3160–M3169 (2013). 3 H. Atsumi, J. Nucl. Mater. 313–316, 543–547 (2003). 4 M. Warrier, R. Schneider, E. Salonen, and K. Nordlund, Nucl. Fusion 47(12), 1656 (2007). 5 A. Ito, Y. Wang, S. Irle, K. Morokuma, and H. Nakamura, J. Nucl. Mater. 390–391, 183–187 (2009). 6 D. C. Elias, R. R. Nair, T. M. G. Mohiuddin, S. V. Morozov, P. Blake, M. P. Halsall, A. C. Ferrari, D. W. Boukhvalov, M. I. Katsnelson, A. K. Geim, and K. S. Novoselov, Science 323(5914), 610–613 (2009). 7 D. W. Boukhvalov and M. I. Katsnelson, Phys. Rev. B 78(8), 085413 (2008). 8 L. Jeloaica and V. Sidis, Chem. Phys. Lett. 300(1–2), 157–162 (1999). 9 L. Hornekær, A. Baurichter, V. V. Petrunin, D. Field, and A. C. Luntz, Science 302(5652), 1943–1946 (2003). 10 X. Sha and B. Jackson, Surf. Sci. 496(3), 318–330 (2002). 11 S. Cazaux and A. G. G. M. Tielens, Astrophys. J. 604(1), 222 (2004). 12 Z. H. Zhu, G. Q. Lu, and F. Y. Wang, J. Phys. Chem. B 109(16), 7923–7927 (2005).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00

134701-9 13 N.

Cortés-Arriagada et al.

Rougeau, D. Teillet-Billy, and V. Sidis, Chem. Phys. Lett. 431(1–3), 135–138 (2006). 14 S. Casolo, O. M. Løvvik, R. Martinazzo, and G. F. Tantardini, J. Chem. Phys. 130(5), 054704 (2009). 15 Y. Wang, H.-J. Qian, K. Morokuma, and S. Irle, J. Phys. Chem. A 116(26), 7154–7160 (2012). 16 M. Bonfanti, R. Martinazzo, G. F. Tantardini, and A. Ponti, J. Phys. Chem. C 111(16), 5825–5829 (2007). 17 G. M. Psofogiannakis and G. E. Froudakis, J. Phys. Chem. C 113(33), 14908–14915 (2009). 18 F. A. Bulat, J. S. Burgess, B. R. Matis, J. W. Baldwin, L. Macaveiu, J. S. Murray, and P. Politzer, J. Phys. Chem. A 116(33), 8644–8652 (2012). 19 J. Murray, Z.-I. Shields, P. Lane, L. Macaveiu, and F. Bulat, J. Mol. Model. 19(7), 2825–2833 (2013). 20 E. Aréou, G. Cartry, J.-M. Layet, and T. Angot, J. Chem. Phys. 134(1), 014701 (2011). 21 Y. Miura, H. Kasai, W. Agerico Diño, H. Nakanishi, and T. Sugimoto, J. Phys. Soc. Jpn. 72(5), 995–997 (2003). 22 R. Martinazzo, S. Casolo, and G. F. Tantardini, in Physics and Applications of Graphene, edited by M. Sergey (InTech, 2011). 23 A. Toro-Labbé, J. Phys. Chem. A 103(22), 4398–4403 (1999). 24 B. Herrera and A. Toro-Labbé, J. Phys. Chem. A 111(26), 5921–5926 (2007). 25 E. Echegaray and A. Toro-Labbé, J. Phys. Chem. A 112(46), 11801–11807 (2008). 26 S. Giri, E. Echegaray, P. W. Ayers, A. S. Nuñez, F. Lund, and A. ToroLabbé, J. Phys. Chem. A 116(40), 10015–10026 (2012). 27 S. Gutiérrez-Oliva, S. Díaz, A. Toro-Labbé, P. Lane, J. S. Murray, and P. Politzer, Mol. Phys. 112(3–4), 349–354 (2013). 28 R. Inostroza-Rivera, B. Herrera, and A. Toro-Labbe, Phys. Chem. Chem. Phys. 16, 14489–14495 (2014). 29 J. Martínez and A. Toro-Labbé, Chem. Phys. Lett. 392(1–3), 132–139 (2004). 30 P. Jaque and A. Toro-Labbé, J. Phys. Chem. A 104(5), 995–1003 (2000). 31 P. Politzer, J. V. Burda, M. C. Concha, P. Lane, and J. S. Murray, J. Phys. Chem. A 110(2), 756–761 (2006). 32 A. Toro-Labbé, S. Gutiérrez-Oliva, J. S. Murray, and P. Politzer, Mol. Phys. 105(19–22), 2619–2625 (2007).

J. Chem. Phys. 141, 134701 (2014) 33 P.

Geerlings, F. De Proft, and W. Langenaeker, Chem. Rev. 103(5), 1793– 1874 (2003). 34 R. Parr and W. Yang, Density Functional Theory of Atoms and Molecules (Oxford University Press, New York, 1989). 35 M. Galván, A. Vela, and J. L. Gazquez, J. Phys. Chem. 92(22), 6470–6474 (1988). 36 S. F. Boys and F. Bernardi, Mol. Phys. 19(4), 553–566 (1970). 37 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77(18), 3865– 3868 (1996). 38 D. Cortés-Arriagada, L. Sanhueza, and K. Wrighton, Int. J. Quantum Chem. 113(15), 1931–1939 (2013). 39 D. Cortés-Arriagada, L. Sanhueza, and M. Santander-Nelli, J. Mol. Model. 19(9), 3569–3580 (2013). 40 S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132(15), 154104 (2010). 41 S. Grimme, S. Ehrlich, and L. Goerigk, J. Comput. Chem. 32(7), 1456– 1465 (2011). 42 E. R. Johnson and A. D. Becke, J. Chem. Phys. 123(2), 174104 (2005). 43 R. Battiti and F. Masulli, in International Neural Network Conference (Springer, Netherlands, 1990), pp. 757–760. 44 F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2(1), 73–78 (2012). 45 M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., Gaussian 09, Revision B.01, Gaussian, Inc., Wallingford, CT, 2009. 46 E. Ghio, L. Mattera, C. Salvo, F. Tommasini, and U. Valbusa, J. Chem. Phys. 73(1), 556–561 (1980). 47 T. Zecho, A. Güttler, X. Sha, B. Jackson, and J. Küppers, J. Chem. Phys. 117(18), 8486–8492 (2002). 48 S. Sakong and P. Kratzer, J. Chem. Phys. 133(5), 054505 (2010). 49 N. Kengo and I. Akira, in Graphene Simulation, edited by J. Gong (InTech, 2011). 50 A. Fukushima, A. Sawairi, K. Doi, M. Senami, L. Chen, H. Cheng, and A. Tachibana, J. Phys. Soc. Jpn. 80(7), 074705 (2011). 51 L. F. Huang, M. Y. Ni, G. R. Zhang, W. H. Zhou, Y. G. Li, X. H. Zheng, and Z. Zeng, J. Chem. Phys. 135(6), 064705 (2011). 52 See supplementary material at http://dx.doi.org/10.1063/1.4896611 for extra tables and figures related to this article.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 137.149.200.5 On: Sat, 22 Nov 2014 12:36:00