The binding of atomic hydrogen on graphene from density functional theory and diffusion Monte Carlo calculations

In this work density functional theory (DFT) and diffusion Monte Carlo (DMC) methods are used to calculate the binding energy of a H atom chemisorbed on the graphene surface. The Perdew-Burke-Ernzerhof (PBE) value of the binding energy is about 20% larger in magnitude than the diffusion Monte Carlo result. The inclusion of exact exchange through the use of the Heyd-Scuseria-Ernzerhof (HSE) functional brings the DFT value of the binding energy closer in line with the DMC result. It is also found that there are significant differences in the charge distributions determined using PBE and DMC approaches.


I. INTRODUCTION
The unique electronic, optical, and transport properties of graphene make it an important system for a wide range of applications, many of which involve or are impacted by the adsorption of atoms or molecules.
To bring these applications to fruition, a deeper understanding of the interaction of atoms and molecules with graphene is required, and, not surprisingly, this has been the subject of several experimental and theoretical studies. [1][2][3][4][5][6][7][8][9][10][11][12][13] The adsorption of H atoms on graphene has been the subject of multiple studies. 3,[10][11][12][13][14] It is known that there is both a weakly absorbed state in which barriers for diffusion are small and a much more strongly bound chemisorbed state. 15,16 Chemisorbed H atoms open up the band gap and allow for tuning of electronic properties. 17 It has been demonstrated that even a single chemisorbed hydrogen atom causes an extended magnetic moment in the graphene sheet. 18,19 On the other hand, there is evidence that given the ready diffusion of H in the physisorbed state, the H atoms tend to pair up on the surface leading to non-magnetic species. 20 Finally, interest in the hydrogen/graphene system has also been motivated by the potential use of graphene and graphitic surfaces for hydrogen storage. 9 The majority of computational studies of adsorption of atoms and molecules on graphene have employed density functional theory (DFT), primarily due to its favorable scaling with system size, allowing for the treatment of larger periodic structures. However, a reliable theoretical description of interactions at the graphene surface has proven to be challenging for DFT. [1][2][3]21 In recent years considerable progress has been made in extending correlated wave function methods to periodic systems. [22][23][24][25][26][27] Among these methods, the diffusion Monte Carlo (DMC) 28 method, which is a real-space stochastic approach to solving the many-body Schrödinger equation is particularly attractive given its low scaling with the number of electrons and high parallelizability. DMC also has the advantages of being systematically improvable and being much less sensitive to the basis set employed than methods that work in the space of Slater determinants. DMC has been used to describe the adsorption of various species on graphene including O 2 4 , a water molecule 5,29 , and a platinum atom. 6 In a study of a physisorbed H atom on graphene, Ma et al. found that different DFT functionals gave binding energies ranging from 5 to 97 meV, while DMC calculations gave a value of only 5 ± 5 meV. 3 Various DFT calculations utilizing the Perdew-Burke-Ernzerhof (PBE) 30 and Perdew-Wang (PW91) 31 functionals predict the chemisorbed H atom species to be bound by 480 to 1,440 meV. [32][33][34][35][36][37][38][39][40] However, this large spread is primarily a result of some calculations employing small supercells resulting in an unphysical description of the low-coverage situation, too small a k-point grid, or small atom-localized basis sets that do not adequately describe the binding and introduce large basis set superposition error (BSSE). In the present work, we use the DMC method to calculate the binding energy of H to graphene in the chemisorbed state.

II. METHODS
All calculations reported in this study used a 5x5x1 supercell of graphene, as it was large enough to make inconsequential the interaction between periodic images of the adsorbed hydrogen atom and to assure that there are essentially unperturbed C atoms between the buckled regions in adjacent images in the x and y directions. The geometries of graphene, both pristine and with a hydro-gen adsorbate, were provided by Kim et al., 41 and were obtained using the PBE+D3 DFT method. 30,42 For all systems, a vacuum spacing of 16Å was used.

A. Density functional theory calculations
The single particle orbitals used in the trial wave functions for variational Monte Carlo (VMC) and DMC calculations were calculated using the PBE functional with the core-correlated electron core potential (ccECP) 43,44 pseudopotentials and a plane wave basis with an energy cutoff of 3,400 eV. Monkhorst-Pack k-point grid meshes 45 were employed with a 13.6 meV Marzari-Vanderbilt-DeVita-Payne cold smearing of the occupations. 46 The PBE results were converged at a 6x6x1 k-point grid.
In addition to the PBE calculations used to generate the trial wave functions for DMC, DFT calculations were carried out with the PBE0 47 and Heyd-Scuseria-Ernzerhof (HSE) functionals 48 to determine if inclusion of exact exchange proves important for the adsorption energy. Due to the inclusion of exact exchange, these calculations would be computationally demanding in a plane wave basis, particularly with the high energetic cutoff required by the ccECP pseudopotential. For this reason, they were carried out all-electron with the POB-TZVP Gaussian type orbital (GTO) basis set. 49 Due to the use of GTOs, these calculations suffer from basis set superposition error (BSSE), for which we use Grimme's geometry-dependent counterpoise correction. 50,51 For the PBE0 and HSE, a 12x12x1 k-point grid was used to assure well converged energies.

B. Quantum Monte Carlo calculations
DMC is a projector quantum Monte Carlo (QMC) method, solving the Schrödinger equation in imaginary time τ = it; any initial state |ψ , that is not orthogonal to the true ground state |φ 0 , will evolve to the ground state in the long time limit. When dealing with Fermionic particles, the DMC method requires the use of the fixed-node approximation 52 to maintain the antisymmetric property of the wave function. For efficient sampling and to reduce statistical fluctuations, we use a Slater-Jastrow trial wave function fixing the nodes through a Slater determinant comprised of single-particle orbitals, which, in this work, are expanded in a B-spline basis. The Jastrow factor is a function that reduces the variance by explicitly describing dynamic correlation. The Jastrow factor contains terms for one-body (electron-ion), two-body (electronelectron) and three-body (electron-electron-ion) interactions. In this study, 10 parameters were employed per spin-channel and the cutoff was fixed to the Wigner-Seitz radius of the simulation cell. The parameters in the Jastrow function were optimized with the linear method 53 using VMC. To reduce the cost of the DMC calculations as well as to reduce the fluctuations near the ionic core regions, ccECP pseudopotentials were used to replace the core electrons. 43,44 The ccECP pseudopotentials were designed to be used with high-accuracy manybody methods such as DMC. The non-local effects due to the ccECP pseudopotentials were addressed using the determinant-localization approximation along with the t-moves method (DLTM). 54,55 Finite size effects were addressed using twist averaging, and symmetry unique twist angles were used. 56 The DMC calculations were performed using the branching scheme proposed by Zen et al. (ZSGMA) 57 with a population control target of 8,192 walkers and a time step of 0.005 a.u., which represented a balance between computational cost and finite timestep error in previous work. 58 We define the binding energy as, where E dgr+H is the energy of the distorted graphene sheet with a chemisorbed atomic hydrogen, E H is the energy of a hydrogen atom, and E gr is the energy of a pristine graphene sheet. In the chemisorbed state, the hydrogen atom bonds directly over a carbon atom, causing this carbon to be pulled out of the sheet towards the hydrogen. 59,60 The adjacent carbons are also pulled in the direction of the hydrogen leading to a distorted graphene sheet. The plane wave DFT calculations were carried out with the QUANTUM ESPRESSO version 6.3 code. [61][62][63] The Gaussian basis DFT calculations were carried out with CRYSTAL17, 64,65 save for the HSE calculation of the lone hydrogen atom which was carried out in NWChem version 6.8 66 using the same basis as the calculations in CRYSTAL17. The QMC calculations were carried out using the QMCPACK code, with the workflow between QUANTUM ESPRESSO and QMCPACK managed by Nexus. [67][68][69] Figures 1 and 3 were rendered with matplotlib 70 and the density plots were generated using VESTA. 71

III. RESULTS & DISCUSSION
A. Binding energy  Table I contains a summary of the binding energies of a hydrogen atom chemisorbed on graphene from this work and selected values from previous publications using the PW91 and PBE functionals. These literature values range from -790 to -z980 meV. However, this wide spread is largely caused by (1) the use in some studies of small supercells for which there are sizable interactions between the CH groups in adjacent cells, and (2) the use in some studies of small atom-centered basis sets without corrections for BSSE. Our calculations with the PBE functional in conjunction with a plane wave basis set give a binding energy of -821 meV. This should be contrasted with our -691 ± 19 meV DMC result. There are several possible sources for the difference between the PBE and DMC values of the binding energy. These include errors in the DFT calculations due to self interaction and planar graphene having more multiconfigurational character than H/graphene, with this being better described with DMC than with PBE. The PBE binding energy is 51 meV lower in magnitude in the plane wave than in the GTO basis set when the same k-point grid is used, and this value is used as a correction for the basis set incompleteness error for the results with other functionals in Table I. The calculations in the GTO basis set give a slightly smaller in magnitude binding energy with PBE0 than with PBE. However, with HSE, we obtain a binding energy 77 meV smaller in magnitude than the PBE result. Applying the correction for the basis set incompleteness error, we obtain -800 meV for the PBE0 binding energy and -743 meV for the HSE binding energy, with the latter being in reasonable agreement with the DMC result of -691 meV.

B. Binding density
It is instructive to examine the change in the electron density associated with the binding of the H atom to the distorted graphene as determined from the PBE and DMC calculations. The density change is given by where ρ H is the charge density of the hydrogen atom, and ρ dgr+H and ρ dgr are the charge densities of the distorted graphene sheet with and without hydrogen, respectively. The ρ b density differences for both DMC and PBE are shown in Figure 2. The dark blue and gold regions represent a loss and gain of electron density, respectively. As expected, there is a shift in electron density from the carbon atom participating in the carbon-hydrogen bond as well as to the three adjacent carbon atoms. These qualitative changes in the density are consistent with previous theoretical and experimental studies. 59,60 The rehybridization from sp 2 to sp 3 of the carbon participating in the CH bond and the weakening of the π bonds due to the distortion of the graphene lead to the electron density shift. The change in the charge distribution is similar for PBE and DMC, with the most noticeable difference being a greater increase of density at remote C atoms in the DMC than in the PBE calculations. Despite the qualitative agreement between the binding density resulting from DMC and PBE, we find that the density differences for the individual systems can offer further insight to the performance of the two methods, as will be seen in the next subsection.

C. Charge density differences between DMC and PBE
In this section, the difference between the DMC and PBE charge densities for distorted graphene with the adsorbed hydrogen atom as well as for pristine planar graphene without the adsorbed hydrogen atom are considered. The charge density difference for each system is calculated according to FIG. 2. Change of the electron density due to the adsorption of the H atom to the distorted graphene sheet (Eq. 2). ρ b from PBE calculations is shown from an oblique angle (A) and aligned along the c axis (B). ρ b from DMC calculations (C) and (D) is shown from the same perspectives. Gold and blue represent a gain and loss of electron density, respectively. Note that there is a region of increased charge density at the C-H bond that is enveloped by a region of loss in the charge density. The binding density was visualized using an isovalue of 2.8×10 −5 for DMC and 3.9×10 −5 for PBE, in both cases capturing 95% of the differential charge density.
where ρ DM C system is the DMC charge density of a given system (either distorted graphene with the adsorbed hydrogen or pristine graphene) and ρ P BE system is the corresponding PBE charge density. ∆ρ gr and ∆ρ dgr+H are reported in Figure 3 along the 110 slice through the unit cell, which captures the carbon-hydrogen bond. From the top-down perspective in Figure 1, the 110 lattice plane bisects the cell diagonally through the longer of the two diagonals. In Figure 3, blue represents areas where the PBE density is larger, while gold areas represent areas where the DMC density is larger. The DMC density, in comparison with the PBE density, has greater weight in the bonding region between atoms. This is the case for both the planar graphene without hydrogen and the system with hydrogen chemisorbed to graphene.
Even though there are significant differences between the PBE and DMC densities for both systems, the difference is similar in the two systems, consistent with it not introducing a large error in the PBE value of the binding energy.

IV. CONCLUSIONS
Calculations of the binding energy of a hydrogen atom on a graphene sheet were carried out using various DFT methods and with DMC. The DMC calculations provide a benchmark value of the binding energy.
Our best estimate of the binding energy from DMC calculations is -691 ± 19 meV. The PBE result obtained with a plane-wave basis set gives a binding energy about 20% larger in magnitude than the DMC result. The global hybrid functional, PBE0, gives a binding energy close to that of PBE. In comparison, HSE, a range-separated hybrid functional, gives a smaller binding energy of -743 meV, after a correction applied for the basis set incompleteness error, and is much closer to the value of the DMC. Interestingly, there are significant differences in the DMC and PBE charge densities of both graphene and H/graphene. Most studies of surface adsorption comparing DMC and DFT results have focused on structures and binding energies. Our work suggests that in terms of understanding performance of DFT methods it is also important to compare the charge densities. The importance of the characterization of electron densities by density functionals was the topic of a recent paper by Medevev et al. 72  FIG. 3. Visualization of the difference of PBE and DMC densities sliced along the 110 lattice plane of the unit cell for the graphene sheet, ∆ρgr, (top) and H adsorbed onto graphene, ∆ρ dgr+H , (bottom). The abscissa represents traversing the 110 plane in fractional coordinates, while the ordinate represents traversing the c axis in fractional coordinates. Blue regions represent places where the PBE density is larger, while the gold color represents regions where the DMC density is larger.

DATA AVAILABILITY STATEMENT
The supplementary material document includes the total energies and error bars for the quantum Monte Carlo calculations, the total energies for the DFT calculations, and details of the convergence of the DFT total energies with respect to the k-point grid and kinetic energy cutoff of the plane wave basis. The data that support the findings of this study are openly available on the Materials Database Facility at https://acdc.alcf.anl.gov/mdf/ detail/dumi_dmc_hgraphene_v1.3, with the following DOI: 10.18126/s1wc-tya. ences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (IN-CITE) program. The DMC and PW-DFT calculations used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC02-06CH11357. The DFT calculation using Gaussian orbitals were carried out on computing resources in the University of Pittsburgh's Center for Research Computing. K.D.J. acknowledges NSF (CBET-2028826) for partial support of this work. S.U. was supported in part by the Pittsburgh Quantum Institute (PQI) Graduate Quantum Leader Award.