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

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–13 

The adsorption of H atoms on graphene has been the subject of multiple studies.^3,10–14 It is known that there are both a weakly absorbed state in which barriers for diffusion are small and a much more strongly bound chemisorbed state,^15,16 which is the focus of this work. Chemisorbed H atoms open up the bandgap and allow for tuning of electronic properties.¹⁷ It has been demonstrated that even a single chemisorbed hydrogen atom causes an extended magnetic moment in a graphene sheet or its derivatives.^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.²⁰ Finally, interest in the hydrogen/graphene system has also been motivated by the potential use of graphene and graphitic surfaces for hydrogen storage.⁹ Despite the interest in H chemisorbed on graphene, we are unaware of an experimental values of the binding energy.

Most 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 large periodic structures. However, a reliable theoretical description of interactions at the graphene surface has proven to be challenging for DFT.^1–3,21 In recent years, considerable progress has been made in extending correlated wave function methods to periodic systems.^22–27 Among these methods, the diffusion Monte Carlo (DMC)²⁸ 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 its energy being much less sensitive to the basis set employed than methods that work in the space of Slater determinants. In DMC calculations, the atomic basis set is important only to the extent that it impacts the nodal surface. DMC has been used to describe the adsorption of various species on graphene including O₂,⁴ a water molecule,^5,29 and a platinum atom.⁶ 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.³ Various DFT calculations utilizing the Perdew–Burke–Ernzerhof (PBE)³⁰ and Perdew–Wang (PW91)³¹ functionals predict the chemisorbed H atom species to be bound by 480 to 1440 meV.^32–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 errors (BSSEs). In the present work, we use the DMC method to calculate the binding energy of H to graphene in the chemisorbed state.

All calculations reported in this study used a 5 × 5 × 1 supercell 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 chemisorbed H atom, were provided by Kim et al.⁴¹ and were obtained using the PBE+D3 DFT method.^30,42 For all systems, a vacuum spacing of 16 Å was used.

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 correlation consistent electron core potential (ccECP)^43,44 pseudopotentials and a plane wave basis with an energy cutoff of 3400 eV. Monkhorst–Pack k-point grid meshes⁴⁵ were employed with a 13.6 meV Marzari–Vanderbilt–DeVita–Payne cold smearing of the occupations.⁴⁶ The PBE results were converged at a 6 × 6 × 1 k-point grid to 1 meV for graphene and graphene with an adsorbed hydrogen atom. The hydrogen atom trial function was generated using a 1 × 1 × 1 k-point grid. Convergence studies can be found in Tables S1 and S2 of the supplementary material.

In addition to the PBE calculations used to generate the trial wave functions for DMC, DFT calculations were carried out with the PBE0⁴⁷ and Heyd–Scuseria–Ernzerhof (HSE) functionals⁴⁸ 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 prohibitive in a plane wave basis, particularly with the high energetic cutoff required by the ccECP pseudopotential. For this reason, they were carried out with the all-electron POB-TZVP Gaussian-type orbital (GTO) basis set.⁴⁹ Due to the use of GTOs, these calculations suffer from basis set superposition errors (BSSEs), which we corrected using Grimme’s geometry-dependent counterpoise correction scheme.^50,51 This correction resulted in a 113 meV reduction in the magnitude of the binding energy when using the PBE0 functional. For the PBE0 and HSE calculations, a 12 × 12 × 1 k-point grid was used to assure binding energies converged to within 2 meV. Convergence data are supplied in Table S3 of the supplementary material.

The plane wave DFT calculations were carried out with the QUANTUM ESPRESSO version 6.3 code.^52–54 The Gaussian basis DFT calculations were carried out with CRYSTAL17,^55,56 save for the HSE calculation of the isolated hydrogen atom that was carried out using NWChem version 6.8⁵⁷ using the same basis as the calculations in CRYSTAL17.

DMC is a projector quantum Monte Carlo (QMC) method, solving the Schrödinger equation in imaginary time τ = it; any initial state |ψ⟩, which is not orthogonal to the true ground state |ϕ₀⟩, 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⁵⁸ 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 (electron–electron), and three-body (electron–electron–ion) interactions. The one- and two-body terms were described with spline functions,⁵⁹ while the three-body terms were represented by polynomials.⁶⁰ Ten parameters were used for the one-body terms per atom type, and ten parameters were employed per spin-channel for the two-body terms. The cutoffs on the one- and two-body terms were fixed to the Wigner–Seitz radius of the simulation cell. The three-body terms were comprised of 26 parameters per term with a cutoff of 10 bohrs. The parameters in the Jastrow factor were separately optimized for each geometry with the linear method⁶¹ 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 many-body methods. The non-local effects due to the pseudopotentials were addressed using the determinant-localization approximation along with the t-moves method (DLTM).^62,63 Finite size effects were addressed using twist averaging.⁶⁴ The twist angles were chosen to be the symmetry unique points of the 6 × 6 × 1 k-point grid shifted by half a grid step away from the gamma point in each direction.

The DMC calculations were performed using the branching scheme proposed by Zen et al. (ZSGMA)⁶⁵ with a population control target of 8192 walkers and a time step of 0.005 a.u., which represented a balance between computational cost and finite timestep error in previous work.⁴ 

We define the binding energy as the following:

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 toward the hydrogen.^66,67 The adjacent carbons are also pulled in the direction of the hydrogen leading to a distorted graphene sheet.

The QMC calculations were carried out using the QMCPACK code, with the workflow between QUANTUM ESPRESSO and QMCPACK managed by Nexus.^69–70 Figures 1 and 3 were rendered with matplotlib,⁷¹ and the density plots were generated using VESTA.⁷² 

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 −980 meV. This widespread of binding energies is 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. We note that the inclusion of the D3 dispersion correction with the PBE functional only changes the magnitude of the binding energy by 30 eV.

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 that obtained in 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. Although the 130 meV difference between the plane-wave PBE and DMC values of the binding energy may appear to be small, this energy difference of that magnitude is consistent with an order of magnitude change in the hydrogen evolution current at room temperature on graphene electrodes.⁴¹ 

In order to better understand the origin of the difference in the PBE and HSE H-atom adsorption energies, we also carried out non-self-consistent calculations using PBE densities to evaluate the HSE energies. These calculations gave a binding energy of only 21 meV smaller in magnitude than that obtained from the self-consistent HSE calculations. This demonstrates that the functional is more important than the density in establishing the binding energy. Detailed information can be found in Table S4 of the supplementary material.

Detailed results of the DMC calculations can be found in Tables S5–S7 in the supplementary material.

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. For the QMC density, the density was accumulated during the VMC and DMC calculations, the mixed estimator bias was found to be insignificant and was, thus, not corrected.

The ρ_b density differences for both DMC and PBE are shown in Fig. 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.^66,67 The rehybridization from sp² to sp³ 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 larger increase of density at remote C atoms in the DMC than in the PBE calculations.

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

where $ρsystemDMC$ is the DMC charge density of a given system (either distorted graphene with the adsorbed hydrogen or pristine graphene) and $ρsystemPBE$ is the corresponding PBE charge density. Δρ_gr and Δρ_dgr+H are reported in Fig. 3 along the 110 slice through the unit cell, which captures the carbon–hydrogen bond. From the top-down perspective in Fig. 1, the 110 lattice plane bisects the cell diagonally through the longer of the two diagonals and is indicated by the solid cyan line. In Fig. 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 larger weight in the bonding region between atoms. We note that the HSE density displays similar differences as the PBE density. (Figure S2) of the supplementary material includes a visualization of the DMC–HSE density difference. 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 the density differences not introducing a large error in the PBE value of the binding energy.

Calculations of the binding energy of a hydrogen atom chemisorbed 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 for the basis set incompleteness error, which is much closer to the value from DMC calculations. Interestingly, there are significant differences in the DMC and PBE charge densities of both graphene and H/graphene.

See the supplementary material for the total energies and error bars for the quantum Monte Carlo calculations, the total energies for the DFT calculations, 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, and a comparison of the density difference of DMC–PBE and DMC–HSE.

We thank Dr. Dan Sorescu for helpful discussions and for sharing the coordinates of his calculations. A.B. and H.S. were supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences 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 (INCITE) program. The DMC and plane wave DFT calculations used the resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357. The DFT calculation using Gaussian orbitals was carried out on computing resources in the University of Pittsburgh’s Center for Research Computing. K.D.J. acknowledges the NSF (Grant No. CBET-2028826) for partial support of this work. S.U. was supported, in part, by the Pittsburgh Quantum Institute (PQI) Graduate Quantum Leader Award.

The authors have no conflicts to disclose.

The data that support the findings of this study are openly available in Materials Database Facility at https://acdc.alcf.anl.gov/mdf/detail/dumi_dmc_hgraphene_v1.3, with the following http://doi.org/10.18126/s1wc-tya.