D

Article pubs.acs.org/JPCC

QM/MD Simulations on Graphene Hydrogenation/Deuteration: CxH/ D Formation Mechanism and Isotope Effect Ying Wang,*,† Hujun Qian,‡ Zhijian Wu,† and Stephan Irle§ †

State Key Laboratory of Rare Earth Resource Utilization, Changchun Institute of Applied Chemistry, Chinese Academy of Sciences, Changchun 130022, P. R. China ‡ State Key Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, Jilin University, Changchun 130023, P. R. China § WPI-Institute of Transformative Bio-Molecules and Department of Chemistry, Graduate School of Science, Nagoya University, Nagoya 464-8602, Japan S Supporting Information *

ABSTRACT: In this work, we performed QM/MD simulations to investigate the hydrogenation process on quasi-free-standing graphene. We found that, at a lower hydrogen incident energy of 0.4 eV, chemisorption occurred in principle in random positions on the graphene surface with little preference for the adsorption site, causing H-frustrated adsorption patterns over time. However, during prolonged exposure to hydrogen, a para-adsorption pattern emerged, producing 25% maximum coverage. In contrast, at higher incident energies of 1.0 and 1.5 eV, the H-frustration pattern emerged at all times, and a higher hydrogen maximum coverage of ∼50% or 60% was reached. This is because the higher incident energies of 1.0 and 1.5 eV can overcome the barrier of 0.63 eV for continuing hydrogen chemisorption on the C4H structure. Based on potential energy curves for C4H−H with H, hydrogenation on para neighbors is a barrierless reaction, favoring the formation of armchair-like C2H structures. The isotope deuterium increased the coverage compared with hydrogen at the same incident energy. Furthermore, the magnetic properties of various hydrogen clusters observed in our QM/MD simulations were studied. The current results show that C4H and graphene hydrogenated with even numbers of H atoms are nonmagnetic, whereas graphene with odd numbers of hydrogen atoms have one to three unpaired electrons, showing magnetism. Hydrogenation opened the band gap by 0.3−4.0 eV at various chemisorption patterns. spintronics.9 For instance, the adsorption of a single hydrogen atom can lead to magnetic moments on the neighboring carbon atoms of 1.0 μB.10 A hydrogen pair with both H atoms adsorbed on two nearest-neighbor carbon atoms of graphene, which corresponds to the most energetically stable configuration, results in a nonmagnetic ground state.9−11 On the contrary, hydrogenation on the surface of graphite or multilayer graphene (in Bernal stacking) can lead to all of the atoms being adsorbed on the same sublattice and prefering a ferromagnetic state for low concentrations.12 Graphone,13,14 a modified graphane obtained by removing hydrogen from one side with C/H = 2:1, has been estimated to have a narrow band gap of 0.5 eV and to exhibit the magnetic properties.15,16 Park and co-workers investigated the stability of multilayer graphone on metal surfaces and the conversion process to sp3 carbon films.17 Buonocore et al. found that the total magnetization of graphone adsorbed on copper and quartz surfaces could be reduced by factors of 4 and 2, respectively.18 Recently, singlesided hydrogenation on quasi-free-standing graphene, which

1. INTRODUCTION Since graphene, a single layer of hexagonally arranged sp2hybridized carbon atoms, was discovered by Novoselov et al.,1 it has attracted a great deal of attention and a number of studies because of its unique physicochemical properties. However, pristine graphene is a nonmagnetic zero-gap semiconductor, which prevents graphene from being used for nanoscale electronic devices with on and off states. Therefore, various efforts have been employed to open the electronic band gap of graphene.2,3 Chemical modification can disorder the structure of graphene with functional groups or create derivatives and is widely applied to open the band gap of graphene. Considering all of the addition species, the greatest possible modification is the addition of hydrogen atoms on graphene to build graphane, which has been predicted theoretically.4 Experimentally, Raman studies revealed that hydrogenation can interrupt the bonding system of graphene through the formation of sp3 carbon− hydrogen bonds and further open the band gap.5 Moreover, the band gap can be tuned in hydrogenated bilayer graphene under a perpendicular electric bias.6 It has also been predicted that partial hydrogenation might induce interesting magnetic properties in graphene7,8 with potential applications in © 2017 American Chemical Society

Received: February 20, 2017 Revised: March 29, 2017 Published: April 6, 2017 8480

DOI: 10.1021/acs.jpcc.7b01662 J. Phys. Chem. C 2017, 121, 8480−8489

Article

The Journal of Physical Chemistry C

and spin-polarization effects play an important role in various chemisorption patterns, SDFTB was employed in our quantum chemical molecular dynamics (QM/MD) simulations. This method was successfully applied in our previous simulations of C4H formation and isotope effect.19,39 Standard C−C DFTB parameters were selected from the mio-0−1 set, which is freely available at http://www.dftb.org. We modified the repulsive potential of the existing mio-0−1 C−H DFTB parameters to closely reflect CCSD(T)/cc-pVTZ(using the G2MS40 approximation) quality potential energy curves of hydrogen chemisorbing on pyrene and coronene.41 The details on the fitting of the new C−H parameters will be published elsewhere. Direct SDFTB trajectories were run by calculating analytical energy gradients on the fly with a velocity Verlet integrator, using 10 au (0.242 fs) as the time step (Δt). Compared with the time step for SWCNT, graphene, or fullerene growth (Δt = 1 fs), this short time step was employed because of the presence of light hydrogen atoms in these systems. Initial velocities were assigned by the prior equilibration of graphene at the target temperature (300 K) for 10 ps. The graphene model system was composed of 32 carbon atoms in a planar unit cell (9.856 × 8.596 × 100 Å3). To study the dependence of coverage and chemisorption patterns on the incident energy, hydrogen or deuterium (H/D) atoms were shot in 0.5-ps intervals with incident energies of 0.4, 1.0, and 1.5 eV in a perpendicular direction from randomly selected positions 3 Å above the graphene surface and with randomly selected spin. The incident energy of 0.4 eV was selected because it was just barely enough to overcome the barrier to hydrogen chemisorption according to previous studies.41,42 The incident energy of 1 eV corresponds to the chemisorption peak from 1 to 100 eV. The incident energy of 1.5 eV was chosen to get a high hydrogen coverage. In total, 500 H/D atoms were shot within the simulation time of 250 ps. To some extent, this designed one-sided shooting simulation is equivalent to hydrogenation on quasi-free-standing graphene (epitaxial graphene on a substrate, in which the back side is passivated by the substrate). The simulations were carried out in the NVT ensemble, in which a constant temperature of 300 K was kept using a Nosé−Hoover chain thermostat. For convergence, an electronic temperature of 300 K was used. At each incident energy, 10 trajectory replicas with different initial conditions were performed. The simulated 50 trajectories are denoted as H0.4a−j, H1.0a−j, H1.5a−j, D0.4a−j, and D1.0a−j, respectively. Band structure calculations were performed using the Vienna ab initio simulation package (VASP).43 The interactions between valence electrons and ion cores were treated by Blöchl’s all-electron-like projector augmented wave (PAW) method. The exchange-correlation functional was the generalized gradient approximation with the Perdew−Burke−Ernzerhof functional (GGA-PBE).44 The wave functions at each kpoint were expanded with a plane-wave basis set and a kinetic energy cutoff of 300 eV. The electron occupancies were determined according to the Fermi scheme with an energy smearing of 0.1 eV. Brillouin zone integration was approximated by a sum over special selected k-points using the Monkhorst−Pack method,45 and they were set to 3 × 3 × 1. Geometries were optimized until the energy was converged to 1.0 × 10−6 eV/atom and the force was converged to 0.01 eV/Å. Spin polarization was considered. A semiempirical DFT-D2 force-field approach46,47 including the van der Waals interaction was employed in our calculations.

was distinguished from traditional double-sided hydrogenation and graphone, was achieved by Haberer et al., and the C4H structure was predicted by both experimental and theoretical studies.19 Because of the interaction between the substrate and graphene, the formation of graphane-type structures was prevented, and the hydrogen uptake was reduced. By varying the hydrogen coverage, a tunable band gap ranging from 1 to 3.5 eV was observed.20 It is known that the chemisorption pattern5,21,22 strongly affects the peculiar electronic and magnetic properties.23 However, much less is known from experiments or simulations about the hydrogenation kinetics on monolayer quasi-free-standing graphene. The formation mechanism of stable phases with certain chemisorption patterns and the electronic properties of these partially hydrogenated graphenes are still not clear. In addition, the chemical adsorption of hydrogen atoms on graphite surfaces is also relevant to a broad range of other fields, such as nuclear fusion, hydrogen storage,24,25 and interstellar chemistry. Therefore, it is urgent to study hydrogenation on quasi-free-standing graphene. In this work, we performed quantum mechanical/molecular dynamics (QM/MD) simulations based on the spin-polarized self-consistent-charge density-functional tight-binding (SDFTB) method. We systematically investigated the kinetics of hydrogenation by considering the nonequilibrium factors in our QM/MD simulations and discussed the coverage dependence on the hydrogen incident energy, as well as the isotope effect. Also, the magnetism and charge distribution were systematically determined for a series of different numbers of hydrogenated clusters, providing a necessary supplement for the magnetic and charge properties of hydrogenated graphene. Theoretical evidence for the existence of stable C4H and armchair-like C2H phases is provided. This knowledge will not only aid in understanding the basic aspects of hydrogenated graphene and other surface phenomena, but will also provide useful information for designing other CxH/CxD species for future experiments.

2. COMPUTATIONAL METHODOLOGY The density functional tight-binding (DFTB) method was the central method employed in the present studies. All DFTB calculations were carried out with the DFTB+ program package developed by Frauenheim and co-workers.26−28 DFTB is an approximate density functional theory method based on the tight-binding approach and utilizes an optimized minimal linear combination of atomic orbitals (LCAO) Slater-type all-valence basis set in combination with a two-center approximation for Hamiltonian matrix elements. Currently, three options are available: (1) In the non-self-consistent-charge (NCC) method, the interactions between atomic charges are neglected. This method is suitable for systems in which the charge interactions are not important, such as systems that include only one type of element. It was successfully used in our previous fullerene formation simulations.29,30 (2) In the self-consistent-charge (SCC) method, in contrast to NCC-DFTB, the interactions between the atomic charge are included, and the energy with atomic charges is determined self-consistently. We have used this method successfully to simulate the growth processes of single-walled carbon nanotubes31−33 and graphene.34−38 (3) In the spin-polarized self-consistent-charge method, on-center spin−spin interactions are also included in the energy, and atomic spin densities are determined self-consistently. This is the most expensive option and is approximately 10 times slower than the NCC option. Because hydrogen has only one electron 8481


Article


Figure 1. (a−c) Final snapshots of trajectories (a) H0.4b/c/d (C32H8 = C4H), H1.0b (C32H12 = C2.7H), and H1.5a (C32H16 = C2H) of the simulated hydrogenation of graphene at 0.4, 1.0, and 1.5 eV for H-shooting followed by 250 ps (500 H) of simulation time. (d,e) Final snapshots of trajectories (d) D0.4 g/j (C32D10 = C3.2D) and (e) D1.0i (C32D19 = C1.7D) of the simulated deuteration of graphene at 0.4 and 1.0 eV for D after 250 ps (500 D) of simulation time. The cyan and red balls are carbon and hydrogen/deuterium atoms, respectively. Black and red circles highlight the hexagons without H/D adsorption (aromatic rings) and H/D back-adsorption, respectively.

Figure 2. Trajectory H0.4b: snapshots of C4H structure formation.

3. RESULTS AND DISCUSSION 3.1. Hydrogen Chemisorption Kinetics. 3.1.1. CxH Formation Mechanism (x ≤ 4). Figure 1 shows the final snapshots (both top view and side view) of five representative trajectories at 0.4, 1.0, and 1.5 eV for H and D addition. The other final snapshots are shown in Figures S1−S5 (Supporting Information). It can be clearly seen that, upon hydrogenation, the attacked carbon atoms are pulled out of the surface, and sp3-hybridized carbon atoms are formed with a rippling structure: The “mountain” carbon atom is connected to a hydrogen atom, and the “valley” carbon atom is not linked to hydrogen, consistent with the previous observation.48 These rippled structures dramatically change the electronic properties of graphene. We also found that all final structures were hydrogen-frustrated, except for the C4H structure, in which all the H atoms were adsorbed on the para position. This H frustration is in agreement with the suggestions of Flores et al.21

In the following sections of this article, the kinetics and formation mechanism of the CxH rippled structures are discussed. First, we considered trajectory H0.4b as representative to demonstrate the formation mechanism of a non-H-frustrated structure (C4H, all H atoms in para positions; see Figure 2), which has been confirmed by many theoretical and experimental observations.19,49 Because the reflection process primarily responded to energy transfer and the original H/D chemisorption patterns were less changed, we show only snapshots of the chemisorption and H2 elimination processes in Figure 2 and Movie S1. It can be seen that the first hydrogenation occurred at 6.5 ps. When the 13th hydrogen atom was shot (its initial position is highlighted with a star), it gradually approached the surrounding carbon atoms and then adsorbed on top of the adjacent carbon atom by pulling it out of the surface (the final position of the H atom is indicated by a blue sphere). Similar behaviors of the chemisorption of the 8482


Article


Figure 3. (a) Schematic of the adsorption sites on C4H−H. (b) Potential energy curves at the SDFTB level for the adsorption of one H atom on pure graphene (C + H), C4H (C4H + H), and C4H−H (C4H−H + H). A, B, C, D, and E correspond to the adsorption sites depicted in panel a, and T and S correspond to the triplet and singlet states, respectively.

Figure 4. (a) Coverage [ratio of H(D) to C] as a function of simulation time at various incident energies for H and D. Averaged from 10 trajectories. (b) Dependence of coverage on the incident energy, where green and red lines/points correspond to average/individual-trajectory coverages for D and H, respectively.

the instability of graphone upon hydrogen-atom migration.51 This C4H formation mechanism fortuitously explains the experimental observation of a 25% maximum H coverage on quasi-free-standing graphene.19 Second, based on the stable C4H structure, a H-frustrated structure was observed at higher coverage. To explain the CxH (x < 4) formation mechanism, the potential energy curves (PECs) of single H atoms interacting with graphene (C + H), C4H (C4H + H), and C4H−H (C4H−H + H−X, where X represents the H adsorption sites denoted by A−E in Figure 3) at the SDFTB level are plotted in Figure 3, along with a schematic of the H adsorption sites. It can be seen that the hydrogenation of C4H has a higher barrier of 0.63 eV (see the curve for C4H + H in Figure 3). Thus, the lower incident energy of 0.4 eV was not high enough to overcome this barrier to obtain coverages higher than 25%. Therefore, only the C4H structure was observed at the incident energy of 0.4 eV. However, when a higher incident energy of 1.0 or 1.5 eV was chosen, it was certainly easy to overcome the barrier of 0.63 eV. Therefore, the chemisorption pattern could cross the C4H structure to form H-frustrated structures and reach a higher maximum coverage of ∼60% (see Figure 4). Furthermore, it is noted that the well depth for further hydrogenation on C4H (reaction of C4H + H) was much shallower (almost the same as the sum of the energies of separated C4H and H) than that for graphene hydrogenation (C + H). This indicates that the complex of C4H−H is less stable than C4H. After the formation of doublet C4H−H, further hydrogenation reactions mainly occur on adsorption site D with a singlet owing to the negligible energy barrier (see the curve for C4H−H + H−D−S in Figure 3). This promises the formation of an armchair-like or boat-like C2H structure rather than a chairlike graphone

other H atoms were observed . With increasing time, more and more H atoms (indicated by red spheres surrounded by blue circles) were randomly and continually attached to graphene without any preferred pattern until 58.5 ps. At this time, the 117th H atom shot at the surface (star) “stole” one preadsorbed H atom (red sphere with blue circle) to form the first H2, as shown in Figure 2 (highlighted by a green ellipse), and the chemisorption pattern was gradually altered. It was noted that H2 can also be formed by two preadsorbed para-hydrogens when a H atom bombards graphene, as indicated by the pink ellipse in Figure S6. According to the previous calculations, this H2 formation reaction has a barrier of 1.4 eV with an exothermic reaction energy of 2.04 eV,50 indicating that the formation of H2 from two para-hydrogens is not easy (thermodynamically favorable while kinetically forbidden) at room temperature. However, in our QM/MD simulations, because of the extremely distorted graphene surface and sudden high-energy bombardment, this type of H2 elimination process was observed. After several chemisorption and H2 elimination steps, the chemisorption pattern was reorganized into a paralike structure. Furthermore, we found that, once the para-like pattern formed, the subsequent incoming hydrogen atoms adsorbed exactly on the para positions, as shown at 135 and 137 ps, until a (full para) C4H structure was built up. In addition, after the formation of the para-like structure, the “unsuitable” (i.e., meta) hydrogens could automatically “swim” to the neighboring para-carbons to form the para-patterned structure, as shown by the pink arrow at 104 ps in Figure 2. This behavior is consistent with the lower barrier of 0.5 eV and the larger exothermic energy of 1.17 eV provided by the conversion of the H-pair pattern (meta to para or meta to ortho)50 and also agrees with the previously reported results on 8483


Article


Figure 5. (a−c) Numbers of (left) ortho, meta, and para chemisorbed pairs of H atoms and aromatic rings and (right) H chemisorption, reflection, and H2 formation as functions of simulation time at incident energies of (a) 0.4, (b) 1.0, and (c) 1.5 eV. (d,e) Same for D at (d) 0.4 and (e) 1 eV. Averaged from 10 trajectories.

structure.15−17,51 This is also different from the double-sided hydrogenation process, in which chair graphane is the most stable structure with the lowest binding energy of −12.23 eV per carbon.52 For the additions at the other sites, such as B, C, and E, unlike for D-site addition, the triplet PECs were found to be more preferable than the singlet PEC (barrier heights of 0.63 vs 0.98 eV for site B, 0.62 vs 0.85 eV for site C, and 0.75 vs

1.09 eV for site E), and they almost overlapped with the curves for C4H hydrogenation. This suggests that it was necessary to overcome an energy barrier for such hydrogenation processes. It is noted that, for site-A adsorption, at the triplet state, the reaction barrier increased to be 0.98 eV, less favorable than adsorption at sites B, C, D, and E. This suggests that the formation of four neighboring H atoms in the same hexagon is 8484


Article


Figure 6. For trajectory H0.4b, (a) numbers of unpaired electrons and adsorbed H atoms as functions of simulation time and (b) spin density plots of partially hydrogenated structures. Red and blue represent alpha and beta spins, respectively; cyan and orange spheres are carbon and hydrogen atoms, respectively. The three values under each structure are the formation time of the structure, the number of H atoms, and the number of unpaired electrons (UE). The isovalue is 0.01 au.

immediately be seen that the H/C and D/C coverages increased linearly with slopes of 0.3 and 0.6, respectively. We hope that these results can potentially offer suggestions for experiments to synthesize various CxH species and expand their wide range of applications in various fields. Figure 5 shows the statistics for the reflection, adsorption, and H2/D2 formation processes, along with the numbers of H pairs (meta, ortho, and para), as functions of simulation time. For the 250-ps shooting simulation (in which a total of 500 hydrogen atoms were shot), it can be clearly seen that the reflection process held the highest ratio of more than 80%. Moreover, we found that the adsorption curve was gradually became parallel to the H2/D2 formation curve (see the insets) with increasing time; that is, the difference between the number of H atoms adsorbing and the number of H atoms desorbing became constant, indicating that the simulation reached a dynamic equilibrium. This was also reflected by the constant numbers of ortho, meta, and para H pairs after 50−150 ps, depending on the incident energy and isotope effects. At lower incident energies, such as 0.4 eV, which can just overcome the reaction barrier for hydrogenation on perfect graphene,41 para pairs dominated throughout the whole process, whereas at higher incident energies, especially at 1.5 eV, which is much higher than the hydrogenation barrier on C4H (0.63 eV), meta pairs were preferred. In contrast, for all investigated systems, regardless of the bombardment species (H or D) and the incident energy, ortho pairs always made a lower contribution than the other two types of H pairs. This is consistent with XPS observation s.39 The final ratios of ortho to meta to para pairs at 0.4, 1.0, and 1.5 eV were 0.10:0.22:0.68, 0.19:0.40:0.41, and 0.23:0.46:0.41, respectively, for H, and those at 0.4 and 1.0 eV were 0.13:0.23:0.64 and 0.27:0.47:0.26, respectively, for D. The

more difficult than the formation of two or three neighbors, that is, generation of a four-legged table-like structure was difficult, consistent with a previous study,53 and few were observed in our current simulations, as shown in Figure 1 and Figures S1−S5. 3.1.2. Incident Energy and Isotope Dependence. Figure 4a depicts the dependence of the coverage on the incident energy and the isotope effect. It can be clearly seen that higher incident energies led to higher coverages up to a coverage of 60%. The coverage increased linearly in the initial stage (less than 50 ps) owing to the large number of sp2 carbons available. However, it is noted that, once the incident energy was too high, the penetration induced back-adsorption (see Figure 1e) or second-layer adsorption. This can finally destroy graphene as graphite peeling.54 According to the expression Ekin = 3/2kBT (where kB is the Boltzmann constant), the incident energies of 0.4, 1.0, and 1.5 eV corresponded to exposing graphene to beams of atomic H obtained by cracking H2 at ∼3000, 7500, and 12000 K, respectively. Thus, to obtain a C4H aromatic structure or a C2H armchair-like structure with 25% or 50% H/ C coverage, a cracking H2 temperature of 3000 or 7500 K, respectively, would be required. Furthermore, it can be seen that the isotope effect was obvious and the coverage of D was always higher than that of H regardless of the incident energy. This resulted from the lower D adsorption barrier and the higher D2 desorption barrier compared with the corresponding values for hydrogenation (considering zero-point-energy corrections), as well as the lower velocity of D to fully relax the bombardment energy because of its heavier mass.39 To clearly see the energy dependence and isotope effect, we fit the coverage curves as a function of incident energy for hydrogenation and deuteration, as shown in Figure 4b. It can 8485


Article


Figure 7. (a−h) Optimized structures of (a) C2H-armchair, (b) C2H-chair, (c) C2H-boat, (d) C2H-zigzag, and (e) C2H-twist and the most stable (f) C2.7H, (g) C3.2H, and (h) C4H. (i) Band structures of C2H-armchair. (j) Band gap and adsorption energy averaged per C−H as a function of the number of H atoms for the selected most stable structures.

density distribution of trajectory H0.4b was determined, as shown in Figure 6b. One can clearly see that the adsorption of a single hydrogen atom enhanced the spin densities of ortho- and para-carbon atoms near the H-attached carbon (as seen at 6.5 ps). This is consistent with the first-principles electronic structure calculations.56 However, for multiple hydrogenation, the spin density on the carbon atoms will be strengthened or canceled by the inducement of multi-H adsorption. For instance, at 16 ps, adsorption of ortho-H1 and -H3 induced a higher spin density on C1 (ortho-carbon). At the same time, meta-H2 lowered the spin density of C1 (meta-carbon) and finally led to no obvious spin density on C1. On the contrary, at 27.5 ps, the three adsorbed hydrogen atoms of ortho-H3 and -H4 and para-H1 led to an increase in the spin density on C2. For the other odd-hydrogenated structures (at 58.5−135 ps), the spin density distributions can be estimated by this regular rule similarly to the spin density distribution, the Mulliken atomic charge distribution were also found to obey this “triangle” charge distribution (as highlighted in blue), in agreement with the previous observations,56 as shown in Figure S7. Also, multibody and mixed-charge interactions will redistribute (strengthen or weaken) the charge on each carbon atom. The spin and Mulliken atomic charge distributions induce the favorable adsorption positions (ortho and para) for subsequent hydrogenation. 3.2.2. Structures, Band Gaps, and Stabilities of CxH (x = 4, 3.2, 2.7, and 2). To investigate the morphologies and properties of partially hydrogenated graphene, the typical structures of CxH (x = 2, 2.7, 3.2, and 4) were optimized, as shown in Figure 7. It can immediately be seen that C4H is uniform with conjugation or aromaticity, where two paracarbon atoms connect to hydrogen to form sp3 carbons and the other four carbon atoms are free of hydrogen and have sp2 hybridization. The stability and ring anisotropy of C4H have

results at the incident energy of 0.4 eV were in good agreement with the experimental hydrogenation observations (0.14:0.24:0.62).19,39 The other results require further experimental confirmation. The number of aromatic rings (ARs, defined as hexagons without hydrogenation) decreased rapidly, especially at higher incident energies. This result agrees well with the trend of increasing coverage (see Figure 2a). The number of ARs finally converged to 4−0, where 4 and 2 indicate C4H para structures and C2H structures, respectively, whereas the other values in this range correspond to other Hfrustrated structures. 3.2. Structures and Properties of Partially Hydrogenated Graphene. 3.2.1. Magnetism and Mulliken Charges of Partially Hydrogenated Graphene. It is worth analyzing the relationship between the number of adsorbedd H atoms and the number of unpaired electron (UEs). Because all trajectories exhibited similar trends, here, we consider only trajectory H0.4b as a representative trajectory. In Figure 6a, one can see that the number of UEs had values of −1, 0, 1, and 3 as the number of H atoms increased continuously from zero to eight. The negative values indicate that the number of beta electrons was larger than the number of alpha electrons. It was also found that, for the adsorption of even numbers of H atoms, such as two, four, six, and eight, the singlet state with zero UEs was always kinetically preferred. However, the adsorption of odd numbers of H atoms created magnetism in the hydrogenated graphene. This is same as the observation of hydrogen adsorption on α-graphyne, in which the adsorption of one hydrogen atom leads to a ferromagnetic configuration and the adsorption of the second hydrogen atom eliminates the magnetization.55 The C4H structure with eight H atoms is a nonmagnetic species, which agrees with the observations of Ranjbar et al.56 To analyze the magnetism of hydrogenated graphene containing odd numbers of hydrogen atoms, the spin 8486


Article

The Journal of Physical Chemistry C been proposed on the basis of many experiments57 and theoretical studies with triangle and “star” structures,49,58 as well as our previous studies.19,39 C 3.2 H, C 2.7 H, and C 2H were formed kinetically by continuous H adsorptions on the ortho/meta positions based on the C4H structure. For C2H, according to a previous study,59 five types of structures, zigzag-like (called C2H-zigzag), armchair-like (called C2H-armchair), boat-like (called C2Hboat), chairlike (called C2H-chair), and twist (called C2Htwist), were investigated. Their stabilities were found to be in the order C2H-armchair > C2H-boat = C2H-twist > C2H-zigzag > C2H-chair. Given that the C2H structures were derived from para-C4H, one-half of hydrogenated C2H structures on substrate were more like C2H-armchair. This armchair structure was similar to the C3H2 structure predicted by Gao et al.60 On the contrary, it was less possible to synthesize C2H-chair, at 17 μB, by this hydrogenation method. The band structure of the most stable C2H (C2H-armchair) structure is shown in Figure 7i. It can be clearly seen that this structure has a direct band gap of 2.6 eV, which is smaller than the band gaps of double-sided-hydrogenated chair and boat structures (3.5 eV).52 From Figure 7j, one can see that hydrogenation opened the band gap of graphene by ∼0.3−4.0 eV. It can also be noted that hydrogenated graphene with odd numbers of H atoms resulted in small band gaps whereas those with even numbers of H atoms led to large band gaps. Furthermore, the average adsorption energy per C−H was also found to depend on the number of H atoms. Specifically, the adsorption of even numbers of H atoms was more stable than the attachment of odd numbers of H atoms. This “zigzagshaped” curve based on the number of H atoms is caused by the different spin states of CxH, that is, odd and even numbers of H atoms correspond to triplet and singlet states, respectively. This is similar to Cx adsorption on the Cu(111) surface.38 Also, we noted that C4H featured the minimum stationary point in the entire hydrogenation process, whereas the adsorption energies of the other structures increased as the odd/even numbers of H atoms increased, leading to a continual decrease in the stability.

formation of an armchair-like C2H structure that was different from the previously reported meta-C2H structure (graphone). The magnetism, band gap, and stability of multihydrogenated quasi-free-standing graphene were found to depend on the number of adsorbed hydrogen atoms. Graphene with even numbers of H atoms had no magnetic properties, larger gaps, and higher stabilities, whereas graphene with odd numbers of H atoms always had one or three unpaired electrons, resulting in magnetic properties, lower band gaps, and lower stability.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b01662. Final snapshots of trajectories H0.4a−j, H1.0a−j, H1.5a−j, D0.4a−j, and D1.0a−j for simulated hydrogenated graphene at 0.4, 1.0, and 1.5 eV after 250 ps (500 H/D) of simulation time (Figure S1−S5); snapshots of para structures formed for H0.4c at 0.4 eV (Figure S6); and Mulliken charge plot for all chemisorption clusters in the MD simulation of trajectory H0.4b at the incident energy of 0.4 eV (Figure S7) (PDF) C4H formation in trajectory H0.4b (Movie S1) (AVI)

■

AUTHOR INFORMATION

Corresponding Author

*Phone: +86-0431-85262801. E-mail: [email protected]. ORCID

Ying Wang: 0000-0002-5437-8741 Hujun Qian: 0000-0001-8149-8776 Stephan Irle: 0000-0003-4995-4991 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (21503210, 21521092, 21673220), Jilin Province Natural Science Foundation (20150101012JC), and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (second phase). Part of the computational time was supported by the Performance Computing Center of Jilin University and Changchun Normal University. Y.W. and S.I. acknowledge helpful discussions with Hiroaki Nakamura and Atsushi Ito of the National Institute of Fusion Science in Toki, Japan.

4. CONCLUSIONS In conclusion, QM/MD simulations were carried out to study the kinetics of hydrogen/deuterium chemisorption process on quasi-free-standing graphene. Three processes, namely, chemical adsorption (forward and back adsorption), reflection, and associative H2/D2 elimination, were observed. Chemisorption first occurred in random positions on graphene, so that Hfrustrated adsorption patterns were created. However, during longer H bombardment, a para-adsorption pattern was observed, as H2/D2 elimination and H/D migration reorganized the chemical pattern obtained at the low incident energy of 0.4 eV. This local aromatic para structure (C4H) was the minimum stationary point in the overall kinetic hydrogenation process, with no magnetism and a band gap of ∼3.5 eV, in agreement with previous experimental and theoretical observations. Once higher incident energies were applied, Hfrustrated structures, especially ortho-carbon structures, with higher maximum coverages of ∼50−60% were observed. This can be attributed to the fact that H atoms with higher energies can easily go through the higher reaction barrier of C4H + H (0.63 eV), producing continual hydrogenation on the para-C4H structure. The further H frustration on C4H−H was found to be a barrierless reaction, promising the kinetically favorable

■

REFERENCES

(1) Novoselov, K. S.; Jiang, D.; Schedin, F.; Booth, T. J.; Khotkevich, V. V.; Morozov, S. V.; Geim, A. K. Two-Dimensional Atomic Crystals. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 10451−10453. (2) Zhou, S. Y.; Gweon, G. H.; Fedorov, A. V.; First, P. N.; de Heer, W. A.; Lee, D. H.; Guinea, F.; Castro Neto, A. H.; Lanzara, A. Substrate-Induced Bandgap Opening in Epitaxial Graphene. Nat. Mater. 2007, 6, 770−775. (3) Han, M. Y.; Ozyilmaz, B.; Zhang, Y.; Kim, P. Energy Band-Gap Engineering of Graphene Nanoribbons. Phys. Rev. Lett. 2007, 98, 206805. (4) Sofo, J. O.; Chaudhari, A. S.; Barber, G. D. Graphane: A TwoDimensional Hydrocarbon. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 153401. (5) Elias, D. C.; Nair, R. R.; Mohiuddin, T. M. G.; Morozov, S. V.; Blake, P.; Halsall, M. P.; Ferrari, A. C.; Boukhvalov, D. W.; Katsnelson, M. I.; Geim, A. K.; Novoselov, K. S. Control of Graphene’s Properties

8487


Article

The Journal of Physical Chemistry C by Reversible Hydrogenation: Evidence for Graphane. Science 2009, 323, 610−613. (6) Samarakoon, D. K.; Wang, X. Q. Tunable Band Gap in Hydrogenated Bilayer Graphene. ACS Nano 2010, 4, 4126−4130. (7) Sun, Y.; Yu, G.; Liu, J.; Shen, X.; Huang, X.; Chen, W. Realizing Diverse Electronic and Magnetic Properties in Hybrid Zigzag BNC Nanoribbons Via Hydrogenation. Phys. Chem. Chem. Phys. 2016, 18, 1326−1340. (8) Guo, X.; Wang, C.; Zhou, Y. Electronic and Magnetic Properties of Hydrogenated Graphene Nanoflakes. Phys. Lett. A 2013, 377, 993− 996. (9) Soriano, D.; Muñoz-Rojas, F.; Fernández-Rossier, J.; Palacios, J. J. Hydrogenated Graphene Nanoribbons for Spintronics. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 165409. (10) Soriano, D.; Leconte, N.; Ordejón, P.; Charlier, J.-C.; Palacios, J.-J.; Roche, S. Magnetoresistance and Magnetic Ordering Fingerprints in Hydrogenated Graphene. Phys. Rev. Lett. 2011, 107, 016602. (11) Ferro, Y.; Teillet-Billy, D.; Rougeau, N.; Sidis, V.; Morisset, S.; Allouche, A. Stability and Magnetism of Hydrogen Dimers on Graphene. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 085417. (12) Moaied, M.; Alvarez, J. V.; Palacios, J. J. Hydrogenation-Induced Ferromagnetism on Graphite Surfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 115441. (13) Zhou, J.; Sun, Q. How to Fabricate a Semihydrogenated Graphene Sheet? A Promising Strategy Explored. Appl. Phys. Lett. 2012, 101, 073114−4. (14) Wu, M.; Burton, J. D.; Tsymbal, E. Y.; Zeng, X. C.; Jena, P. Hydroxyl-Decorated Graphene Systems as Candidates for Organic Metal-Free Ferroelectrics, Multiferroics, and High-Performance Proton Battery Cathode Materials. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 081406. (15) Zhou, J.; Wang, Q.; Sun, Q.; Chen, X. S.; Kawazoe, Y.; Jena, P. Ferromagnetism in Semihydrogenated Graphene Sheet. Nano Lett. 2009, 9, 3867−3870. (16) Rudenko, A. N.; Keil, F. J.; Katsnelson, M. I.; Lichtenstein, A. I. Exchange Interactions and Frustrated Magnetism in Single-Side Hydrogenated and Fluorinated Graphene. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 081405. (17) Odkhuu, D.; Shin, D.; Ruoff, R. S.; Park, N. Conversion of Multilayer Graphene into Continuous Ultrathin sp3-Bonded Carbon Films on Metal Surfaces. Sci. Rep. 2013, 3, 3276. (18) Buonocore, F.; Mosca Conte, A.; Lisi, N. Effects of the Substrate on Graphone Magnetism: A Density Functional Theory Study. Phys. E 2016, 78, 65−72. (19) Haberer, D.; Giusca, C. E.; Wang, Y.; Sachdev, H.; Fedorov, A. V.; Farjam, M.; Jafari, S. A.; Vyalikh, D. V.; Usachov, D.; Liu, X.; Treske, U.; Grobosch, M.; Vilkov, O.; Adamchuk, V. K.; Irle, S.; Silva, S. R. P.; Knupfer, M.; Büchner, B.; Grüneis, A. Evidence for a New Two-Dimensional C4H-Type Polymer Based on Hydrogenated Graphene. Adv. Mater. 2011, 23, 4497−4503. (20) Haberer, D.; Vyalikh, D. V.; Taioli, S.; Dora, B.; Farjam, M.; Fink, J.; Marchenko, D.; Pichler, T.; Ziegler, K.; Simonucci, S.; Dresselhaus, M. S.; Knupfer, M.; Büchner, B.; Grüneis, A. Tunable Band Gap in Hydrogenated Quasi-Free-Standing Graphene. Nano Lett. 2010, 10, 3360−3366. (21) Flores, M. Z. S.; Autreto, P. A. S.; Legoas, S. B.; Galvao, D. S. Graphene to Graphane: A Theoretical Study. Nanotechnology 2009, 20, 465704−6. (22) Ryu, S.; Han, M. Y.; Maultzsch, J.; Heinz, T. F.; Kim, P.; Steigerwald, M. L.; Brus, L. E. Reversible Basal Plane Hydrogenation of Graphene. Nano Lett. 2008, 8, 4597−4602. (23) Boukhvalov, D. W.; Katsnelson, M. I.; Lichtenstein, A. I. Hydrogen on Graphene: Electronic Structure, Total Energy, Structural Distortions and Magnetism from First-Principles Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 035427. (24) Meregalli, V.; Parrinello, M. Review of Theoretical Calculations of Hydrogen Storage in Carbon-Based Materials. Appl. Phys. A: Mater. Sci. Process. 2001, 72, 143−146.

(25) Schlapbach, L.; Zuttel, A. Hydrogen-Storage Materials for Mobile Applications. Nature 2001, 414, 353−358. (26) Porezag, D.; Frauenheim, T.; Kohler, T.; Seifert, G.; Kaschner, R. Construction of Tight-Binding-Like Potentials on the Basis of Density-Functional Theory: Application to Carbon. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 51, 12947−12956. (27) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260−7268. (28) Kohler, C.; Seifert, G.; Gerstmann, U.; Elstner, M.; Overhof, H.; Frauenheim, T. Approximate Density-Functional Calculations of Spin Densities in Large Molecular Systems and Complex Solids. Phys. Chem. Chem. Phys. 2001, 3, 5109−5114. (29) Irle, S.; Zheng, G.; Elstner, M.; Morokuma, K. From C2 Molecules to Self-Assembled Fullerenes in Quantum Chemical Molecular Dynamics. Nano Lett. 2003, 3, 1657−1664. (30) Qian, H. J.; Wang, Y.; Morokuma, K. Quantum Mechanical Simulation Reveals the Role of Cold Helium Atoms and the Coexistence of Bottom-up and Top-Down Formation Mechanisms of Buckminsterfullerene from Carbon Vapor. Carbon 2017, 114, 635− 641. (31) Page, A. J.; Yamane, H.; Ohta, Y.; Irle, S.; Morokuma, K. QM/ MD Simulation of SWNT Nucleation on Transition-Metal Carbide Nanoparticles. J. Am. Chem. Soc. 2010, 132, 15699−15707. (32) Wang, Y.; Gao, X.; Qian, H.-J.; Ohta, Y.; Wu, X.; Eres, G.; Morokuma, K.; Irle, S. Quantum Chemical Simulations Reveal Acetylene-Based Growth Mechanisms in the Chemical Vapor Deposition Synthesis of Carbon Nanotubes. Carbon 2014, 72, 22−37. (33) Wang, Y.; Jiao, M.; Wu, Z.; Irle, S. Theoretical Studies on Ethanol Dissociation on Iron Nano-Particles in the Early Stage of Swcnt Growth. J. Phys. Chem. C 2017, 121, 2276−2284. (34) Page, A. J.; Wang, Y.; Li, H.-B.; Irle, S.; Morokuma, K. Nucleation of Graphene Precursors on Transition Metal Surfaces: Insights from Theoretical Simulations. J. Phys. Chem. C 2013, 117, 14858−14864. (35) Wang, Y.; Page, A. J.; Nishimoto, Y.; Qian, H. J.; Morokuma, K.; Irle, S. Template Effect in the Competition between Haeckelite and Graphene Growth on Ni(111): Quantum Chemical Molecular Dynamics Simulations. J. Am. Chem. Soc. 2011, 133, 18837−18842. (36) Jiao, M.; Li, K.; Wang, Y.; Wu, Z. J. Quantum Chemical Molecular Dynamics Studies of Bilayer Graphene Growth on a Ni(111) surface. J. Phys. Chem. C 2015, 119, 12643−12650. (37) Wang, Y.; Page, A. J.; Li, H. B.; Qian, H. J.; Jiao, M. G.; Wu, Z. J.; Morokuma, K.; Irle, S. Step-Edge Self-Assembly During Graphene Nucleation on a Nickel Surface: QM/MD Simulations. Nanoscale 2014, 6, 140−144. (38) Page, A. J.; Mitchell, I.; Li, H. B.; Wang, Y.; Jiao, M.; Irle, S.; Morokuma, K. Spanning the “Parameter Space” of Chemical Vapor Deposition Graphene Growth with Quantum Chemical Simulations. J. Phys. Chem. C 2016, 120, 13851−13864. (39) Paris, A.; Verbitskiy, N.; Nefedov, A.; Wang, Y.; Fedorov, A.; Haberer, D.; Oehzelt, M.; Petaccia, L.; Usachov, D.; Vyalikh, D.; Sachdev, H.; Wöll, C.; Knupfer, M.; Büchner, B.; Calliari, L.; Yashina, L.; Irle, S.; Grüneis, A. Kinetic Isotope Effect in the Hydrogenation and Deuteration of Graphene. Adv. Funct. Mater. 2013, 23, 1628− 1635. (40) Froese, R. D. J.; Humbel, S. p.; Svensson, M.; Morokuma, K. Imomo, a New High-Level G2-Like Method for Large Molecules and ́ Its Applications to Diels, ialder Reactions. J. Phys. Chem. A 1997, 101, 227−233. (41) Wang, Y.; Qian, H. J.; Morokuma, K.; Irle, S. Coupled Cluster and Density Functional Theory Calculations of Atomic Hydrogen Chemisorption on Pyrene and Coronene as Model Systems for Graphene Hydrogenation. J. Phys. Chem. A 2012, 116, 7154−7160. (42) Kerwin, J.; Jackson, B. The Sticking of H and D Atoms on a Graphite (0001) Surface: The Effects of Coverage and Energy Dissipation. J. Chem. Phys. 2008, 128, 084702−084707. 8488


Article

The Journal of Physical Chemistry C (43) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (44) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (45) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188−5192. (46) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104−19. (47) Grimme, S. Semiempirical Gga-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (48) Boukhvalov, D. W.; Katsnelson, M. I. Enhancement of Chemical Activity in Corrugated Graphene. J. Phys. Chem. C 2009, 113, 14176− 14178. (49) Ferro, Y.; Morisset, S.; Allouche, A. Evidence of Hydrogenated Hexamers on Graphite. Chem. Phys. Lett. 2009, 478, 42−44. (50) Hornekaer, L.; Šljivančanin, Ž .; Xu, W.; Otero, R.; Rauls, E.; Stensgaard, I.; Laegsgaard, E.; Hammer, B.; Besenbacher, F. Metastable Structures and Recombination Pathways for Atomic Hydrogen on the Graphite (0001) Surface. Phys. Rev. Lett. 2006, 96, 156104. (51) Boukhvalov, D. W. Stable Antiferromagnetic Graphone. Phys. E 2010, 43, 199−201. (52) Samarakoon, D. K.; Wang, X.-Q. Chair and Twist-Boat Membranes in Hydrogenated Graphene. ACS Nano 2009, 3, 4017− 4022. (53) AlZahrani, A. Z.; Srivastava, G. P. Structural and Electronic Properties of H-Passivated Graphene. Appl. Surf. Sci. 2010, 256, 5783− 5788. (54) Ito, A.; Wang, Y.; Irle, S.; Morokuma, K.; Nakamura, H. Molecular Dynamics Simulation of Hydrogen Atom Sputtering on the Surface of Graphite with Defect and Edge. J. Nucl. Mater. 2009, 390391, 183−187. (55) Drogar, J.; Roknabadi, M. R.; Behdani, M.; Modarresi, M.; Kari, A. Hydrogen Adsorption on the A-Graphyne Using Ab Initio Calculations. Superlattices Microstruct. 2014, 75, 340−346. (56) Ranjbar, A.; Bahramy, M. S.; Khazaei, M.; Mizuseki, H.; Kawazoe, Y. First-Principles Study of Structural Stability, Magnetism, and Hyperfine Coupling in Hydrogen Clusters Adsorbed on Graphene. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 165446. (57) Hornekær, L.; Xu, W.; Otero, R.; Lægsgaard, E.; Besenbacher, F. Long Range Orientation of Meta-Stable Atomic Hydrogen Adsorbate Clusters on the Graphite(0001) Surface. Chem. Phys. Lett. 2007, 446, 237−242. (58) Khazaei, M.; Bahramy, M. S.; Ranjbar, A.; Mizuseki, H.; Kawazoe, Y. Geometrical Indications of Adsorbed Hydrogen Atoms on Graphite Producing Star and Ellipsoidal Like Features in Scanning Tunneling Microscopy Images: Ab Initio Study. Carbon 2009, 47, 3306−3312. (59) Zhou, C.; Chen, S.; Lou, J.; Wang, J.; Yang, Q.; Liu, C.; Huang, D.; Zhu, T. Graphene’s Cousin: The Present and Future of Graphane. Nanoscale Res. Lett. 2014, 9, 26. (60) Gao, B.; Shao, X.; Lv, J.; Wang, Y.; Ma, Y. Structure Prediction of Atoms Adsorbed on Two-Dimensional Layer Materials: Method and Applications. J. Phys. Chem. C 2015, 119, 20111−20118.

8489


D

Recommend Documents