key: cord-0907177-ogq5sihs authors: Dumi, Amanda; Upadhyay, Shiv; Bernasconi, Leonardo; Shin, Hyeondeok; Benali, Anouar; Jordan, Kenneth D. title: The binding of atomic hydrogen on graphene from density functional theory and diffusion Monte Carlo calculations date: 2022-02-01 journal: J Chem Phys DOI: 10.1063/5.0085982 sha: 6de7c323741f0baeb75a8c430b536010280cd074 doc_id: 907177 cord_uid: ogq5sihs 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. 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. 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. 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. 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 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. 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. 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. 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. 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. Computational Investigations of Dispersion Interactions between Small Molecules and Graphene-like Flakes Permeation of chemisorbed hydrogen through graphene: A flipping mechanism elucidated Binding of hydrogen on benzene, coronene, and graphene from quantum Monte Carlo calculations Diffusion Monte Carlo study of O2 adsorption on single layer graphene Adsorption and diffusion of water on graphene from first principles Adsorption of a single pt atom on graphene: spin crossing between physisorbed triplet and chemisorbed singlet states Graphene: Status and Prospects Two-dimensional gas of massless Dirac fermions in graphene Graphene and Graphene-Like Materials for Hydrogen Energy Hydrogen storage: Materials, methods and perspectives Adsorption-Based Hydrogen Storage in Activated Carbons and Model Carbon Structures High-capacity hydrogen storage by metallized graphene Pillared Graphene: A New 3-D Network Nanostructure for Enhanced Hydrogen Storage On the performance of van der Waals corrected-density functional theory in describing the atomic hydrogen physisorption on graphite First-principles study of the structural and energetic properties of H atoms on a graphite (0001) surface DFT investigation of the adsorption of atomic hydrogen on a cluster-model graphite surface Band gap opening in graphene: a short theoretical study Atomic-scale control of graphene magnetism by using hydrogen atoms Hydrogen physisorption channel on graphene: a highway for atomic H diffusion Hydrogen physisorption channel on graphene: a highway for atomic H diffusion Low rank factorization of the Coulomb integrals for periodic coupled cluster theory An optimized twist angle to find the twist-averaged correlation energy applied to the uniform electron gas Dynamical correlation energy of metals in large basis sets from downfolding and composite approaches Towards an exact description of electronic wavefunctions in real solids Local embedding of coupled cluster theory into the random phase approximation using plane waves Toward a systematic improvement of the fixed-node approximation in diffusion Monte Carlo for solids-A case study in diamond Quartic scaling mp2 for solids: A highly parallelized algorithm in the plane wave basis Quantum monte carlo simulations of solids Physisorption of water on graphene: Subchemical accuracy from many-body electronic structure methods Generalized Gradient Approximation Made Simple Unified theory of exchange and correlation beyond the local density approximation Hydrogen pairing on graphene Extended atomic hydrogen dimer configurations on the graphite(0001) surface Irradiation-induced magnetism in graphite: A density functional study Hydrogen storage by spillover on graphene as a phase nucleation process Understanding adsorption of hydrogen atoms on graphene Hydrogen adsorption on boron doped graphene: an ab initio study DFT calculation for adatom adsorption on graphene sheet as a prototype of carbon nanotube functionalization Concentration dependent magnetism induced by hydrogen adsorption on graphene and single walled carbon nanotubes Hydrogen on graphene: Electronic structure, total energy, structural distortions and magnetism from first-principles calculations Real Time Modulation of Hydrogen Evolution Activity of Graphene Electrodes Using Mechanical Strain A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu A new generation of effective core potentials from correlated calculations: 3d transition metal series A new generation of effective core potentials for correlated calculations Special points for Brillouinzone integrations Thermal Contraction and Disordering of the Al(110) Surface Toward reliable density functional methods without adjustable parameters: The PBE0 model Influence of the exchange screening parameter on the performance of screened hybrid functionals Consistent gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations A geometrical correction for the interand intra-molecular basis set superposition error in Hartree-Fock and density functional theory calculations for large systems Geometrical Correction for the Inter-and Intramolecular Basis Set Superposition Error in Periodic Density Functional Theory Calculations Quantum chemistry by random walk: Higher accuracy Alleviation of the Fermion-Sign Problem by Optimization of Many-Body Wave Functions A new scheme for fixed node diffusion quantum Monte Carlo with pseudopotentials: Improving reproducibility and reducing the trialwave-function bias Size-consistent variational approaches to nonlocal pseudopotentials: Standard and lattice regularized diffusion Monte Carlo methods revisited Twist-averaged boundary conditions in continuum quantum monte carlo algorithms Boosting the accuracy and speed of quantum Monte Carlo: Size consistency and time step Diffusion Monte Carlo study of O 2 adsorption on single layer graphene The mechanism of chemisorption of hydrogen atom on graphene: Insights from the reaction force and reaction electronic flux Imaging covalent bond formation by H atom scattering from graphene Quantum ESPRESSO toward the exascale QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials Advanced capabilities for materials modelling with QUANTUM ESPRESSO Quantum-mechanical condensed matter simulations with CRYSTAL and beyond, a long story QMCPACK: an open source ab initio quantum Monte Carlo package for the electronic structure of atoms, molecules and solids QMCPACK: Advances in the development, efficiency, and application of auxiliary field and real-space variational and diffusion quantum Monte Carlo Nexus: A modular workflow management system for quantum simulation codes Matplotlib: A 2d graphics environment VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data Density functional theory is straying from the path toward the exact functional We thank Dr. Dan Sorescu for helpful discussion 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 Sci-