Unravelling H$_2$ chemisorption and physisorption on metal decorated graphene using quantum Monte Carlo

Molecular hydrogen is at the core of hydrogen energy applications and has the potential to significantly reduce the use of carbon dioxide emitting energy processes. However, hydrogen gas storage is a major bottleneck for its large-scale use as current storage methods are energy intensive. Among different storage methods, physisorbing molecular hydrogen at ambient pressure and temperatures is a promising alternative - particularly thanks to tuneable lightweight nanomaterials and high throughput screening methods. Nonetheless, understanding hydrogen adsorption in well-defined nanomaterials remains experimentally challenging and reference information is scarce despite the proliferation of works predicting hydrogen adsorption. In this work, we focus on Li, Na, Ca, and K, decorated graphene sheets as substrates for hydrogen adsorption and compute the most accurate adsorption energies available to date using quantum diffusion Monte Carlo (DMC). Building on our previous insights at the density functional theory (DFT) level, we find that a weak covalent chemisorption of molecular hydrogen, known as Kubas binding, is feasible on Ca decorated graphene according to DMC, in agreement with DFT. This finding is in contrast to previous DMC predictions of the 4H$_2$/Ca$^+$ gas cluster where chemisorption is not favoured. However, we find that the adsorption energy of hydrogen on metal decorated graphene according to a widely-used DFT method is not fully consistent with DMC and the discrepancies are not systematic. The reference adsorption energies reported herein can be used to find better work-horse methods for application in large-scale modelling of hydrogen adsorption. Furthermore, the implications of this work affect strategies for finding suitable hydrogen storage materials and high-throughput methods.

Molecular hydrogen is at the core of hydrogen energy applications and has the potential to significantly reduce the use of carbon dioxide emitting energy processes. However, hydrogen gas storage is a major bottleneck for its large-scale use as current storage methods are energy intensive. Among different storage methods, physisorbing molecular hydrogen in nanomaterials at ambient pressure and temperatures is a promising alternative -particularly in light of the advancements in tuneable lightweight nanomaterials and high throughput screening methods. Nonetheless, understanding hydrogen adsorption in well-defined nanomaterials remains experimentally challenging and therefore, reference information is scarce despite the proliferation of works predicting hydrogen adsorption. In this work, we focus on Li, Na, Ca, and K, decorated graphene sheets as substrates for molecular hydrogen adsorption and compute the most accurate adsorption energies available to date using quantum diffusion Monte Carlo (DMC). Building on our previous insights at the density functional theory (DFT) level, we find that a weak covalent chemisorption of molecular hydrogen, known as Kubas interaction, is thermodynamically feasible on Ca decorated graphene according to DMC, in agreement with DFT. This finding is in contrast to previous DMC predictions of the 4H 2 /Ca + gas cluster (without graphene) where chemisorption is not favoured. However, we find that the adsorption energy of hydrogen on metal decorated graphene according to a typical widely-used DFT method is not fully consistent with DMC and the discrepancies are not systematic. The reference adsorption energies reported herein can be used to find better performing work-horse methods for application in large-scale modelling of hydrogen adsorption. In addition, this work demonstrates the importance of understanding the interaction mechanisms governing adsorption as well as the corresponding absolute adsorption energies. The implications of this work affect strategies for finding suitable hydrogen storage materials and high-throughput methods.

I. INTRODUCTION
The smallest molecule, hydrogen (H 2 ), has the potential to sharply improve sustainable energy production, by replacing fossil fuels in a variety of applications such as vehicular fuel and home heating. For example, water is the only by-product of converting H 2 using polymer electrolyte membrane (PEM) fuel cells and therefore its use would drastically reduce pollutants from vehicles. However, the efficient storage of hydrogen molecules is an outstanding challenge. The most used storage method currently is to pressurize H 2 at high (∼ 700 bar) pressure inside carbon fiber tanks. 1 This simple but expensive route affects the fuel economy of vehicles, detracting considerably from the benefits of hydrogen technology. Alternative storage strategies includes strong (dissociative) chemisorption, also know as the spillover effect, and non-dissociative physisorption, whereby H 2 is adsorbed on a surface as intact gas molecules. The spillover effect is typically associated with high energy barriers for hydrogen release, which also reduces the energy yield from hydrogen. Meanwhile, the spontaneous adsorption of molecular hydrogen within an energy window of −200 to −400 meV per H 2 molecule in a lightweight material, circumvents the need for high pressures or extreme temperatures. 2,3 Our focus is to find materials and interaction mechanisms that bring the H 2 adsorption energy into that range on a lightweight material.
Graphene, carbon nanotubes, and analogous lowdimensional materials have the advantage of being lightweight with high surface areas for adsorption and being made of earth abundant elements. However, experimentally derived adsorption energies and theoretical predictions for H 2 on graphene and carbon nanotubes appear to be very weak, typically at less than −50 meV per H 2 molecule. [3][4][5][6][7][8][9][10][11] , In our previous work, 2 a density functional theory (DFT) based method predicts that metal decorating atoms on graphene, such as Li and Ca, can boost the adsorption energy by over 100 meV -bringing it almost into the useful binding energy window for storage. However, it is well known that DFT methods suffer from the electron delocalization problem and can also be inaccurate for predicting long-range dispersion based interactions. In particular, it has previously been shown that DFT methods wrongly stabilize a covalent bond between hydrogen and group 2 metals with available 3d-states such as Ca, known as Kubas-type binding, 12 in the gas phase. [13][14][15][16] However, it is not known if this inaccuracy carries over to the solid-state system where the metal atoms are supported by graphene. Since Kubas-type binding is predicted to be a very tunable form of chemisorption interaction in low-dimensional materials at the DFT level, we aim to establish whether it is a feasible interaction using a higher-accuracy method in this work. We use a wavefunction based method that has been shown to have benchmark accuracy, namely, quantum diffusion Monte Carlo (DMC). We predict the most promising materials indicated in our previous work with DMC and report the most accurate adsorption energies available to date for H 2 on pristine and metal decorated graphene. In doing so, we establish the propensity of the materials considered here to bind H 2 . We also find that a widely-used DFT method slightly overestimates the binding for alkali metal (Li, Na, and K) decorated graphene and severely overestimates binding on Ca decorated graphene.

II. METHODS AND COMPUTATIONAL SETUP
We established PBE+D3 binding geometries for hydrogen on pristine graphene in Ref. 5 and on Li, Na, Ca, and K decorated graphene in a recent work. 2 PBE is a widely used generalized gradient approximation in DFT 17 and a semi-empirical two and three-body dispersion correction, D3 is added using the zero-damping function. 18 We use the optimized geometries with the most favourable binding energies per H 2 molecule according to PBE+D3: 3H 2 +Li@Gr, 3H 2 +Na@Gr, 4H 2 +K@Gr, and 4H 2 +Ca@Gr, as can be seen in Fig. 1. The adsorption energy is defined as: where E tot nH 2 +M@Gr is the total energy of n number of hydrogen molecules adsorbed on metal decorated graphene (M@Gr) where the metal can be Li, Na, K, or Ca. Correspondingly, E tot M@Gr is the total energy of the metal decorated graphene substrate and E tot H 2 is the total energy of a gas phase hydrogen molecule. Note that the graphene sheet is a (5 × 5) unit cell of graphene and PAW potentials (available for the VASP 5.4.4 software package) with explicit semi-core electrons were used for the metal atoms. [19][20][21][22] In addition to the geometries previously established, we used PBE+D3 to find a physisorption minimum for 4H 2 molecules on Ca@Gr. This structure is also available in the Supplemental Material (SM) and the corresponding H 2 -Ca separation distances are ca. 3.3-3.6 Å with an adsorption energy of −81 meV per H 2 molecule. In order to understand the impact of electron localization at the DFT level, we use the simplified approach to the Hubbard U method, introduced by Dudarev et al., 23 as implemented in VASP v.5.4.4. These heuristic calculations used U values for 3d-states in Ca between 1 to 6 eV, which is a range that is typically applied to metal atoms. Note that we set the J value to zero in all calculations. The DFT-VASP calculations of the Ca@Gr based system undertaken in this work are computed using a 3 × 3 × 1 k-point grid, shifted from Γ by (1/2, 1/2, 0) and 400 eV plane-wave cut-off. This yields only a 2 meV difference in the (Kubas-bound) 4H 2 +Ca@Gr PBE+D3 adsorption energy with respect to our previous work. 2 Previously, we computed the adsorption energy of a single H 2 molecule on pristine graphene 5 using a (3 × 3) unit cell of graphene and two k-points. Fixed-nodes were used with Slater-Jastrow type trial wavefunctions, where the orbitals for the Slater determinant were computed using PWSCF in Quantum Espresso 24,25 . The result was −24 ± 11 meV. Here, we FIG. 1. The PBE+D3 geometries of (a) H 2 on pristine graphene (H 2 +Gr), (b) 3H 2 +Li@Gr, (c) 3H 2 Na@Gr, (d) 4H 2 +K@Gr, and (e) 4H 2 +Ca@Gr. The unit cell is indicated by a blue bounding box. follow the same computational setup but we invoke recent developments, mainly the determinant localization approximation (DLA) in QMC, 26 which avoids a bias in the QMC energy from the use of different Jastrow factors in the limit of small time-steps. We used CASINO v.2.13 27 (with time-steps of 0.03 and 0.05 a.u.) as well as the GPU enabled DMC algorithm in QMCPACK v.13.5.9 28 (with time-steps of 0.015 and 0.05 a.u.) to obtain the hydrogen adsorption energy with stochastic errors of 2-5 meV. The results from the two QMC codes and different time-steps are all fully consistent within the stochastic error bars. See the SM for full details.
The geometries from PBE+D3 (available in the SM) were used to perform DFT calculations with Quantum Espresso v.6.8 24,25 . The unit cells can be seen in Fig. 1. More specifically, input orbitals were computed using the local density approximation (LDA) in Quantum Espresso, with ccECP pseudopotentials 29-32 , a 400 Ry planewave energy cut-off, and a vacuum of ca. 15 Å in the direction perpendicular to graphene. The ccECP pseudopotentials treated 3 electrons explicitly for Li, 9 electrons for Na and K, and 10 electrons for Ca. We tested the effect of using softer-core ccECP pseudopotential for Na and Li at the DFT level. We found that the PBE+D3 adsorption energy per H 2 molecules deviates by 5 meV for 3H 2 +Li@Gr and 7 meV for 3H 2 +Na@Gr due to the softer pseudopotentials.
The orbitals obtained from DFT are used to define the nodal surface of each system in DMC and it is typically insensitive to the choice of DFT functional that is employed, however we also validate it here using two DFT functionals. LDA is a simple and efficient method that can be used for metallic systems and to validate its use here, we computed the LDA and hybrid B3LYP 33 orbitals for a small gas-phase Ca + +4H 2 cluster for comparison. The resulting DMC total energies and interaction energies based on LDA and B3LYP orbitals are statistically indistinguishable (see Appendix Section B for further information). All calculations were spin-polarized. For the Li, Na, and K decorated graphene systems, a k-point mesh of 3 × 3 × 1 centred on the Γ-point was used which is sufficient to converge hydrogen adsorption energies at the DFT level. Using the resulting charge density, a separate non selfconsistent field (NSCF) calculation was performed at a single k-point to produce orbitals for the following QMC calculations. The adsorption energy based on a single k-point is in agreement with the adsorption energy from the full kpoint grid for 3H 2 +Li@Gr, 3H 2 +Na@Gr, and 4H 2 +K@Gr, as can be seen in the SM. However, the adsorption energy in the 4H 2 +Ca@Gr system is more dependent on the k-point grid and for this system we undertook a careful analysis of the occupations at different k-points. We find that a 3 × 3 × 1 kpoint grid shifted from Γ by (1/2, 1/2, 0) in reciprocal units yields a converged adsorption energy at the LDA level, with respect to a 15 × 15 × 1 k-point mesh. Moreover, partial occupations around the Fermi energy are avoided using the shifted 3 × 3 × 1 k-point grid which is preferential as it allows us to perform QMC calculations at separate k-points using only integer occupations. LDA orbitals were obtained using NSCF calculations at each k-point with occupations set according to the those in the LDA SCF calculation in the full k-point grid.
Further information and justification of this protocol based on the electronic structure for Ca@Gr based systems can be found in Appendix Section A and the SM. The planewave orbitals were localized using B-splines 34 and a meshfactor of 0.5 in QMCPACK.
QMC calculations were performed using GPU enabled complex version of QMCPACK v.3.15.9. 28 We optimized one, two, and three-body parameters for the Jastrow factor using variational Monte Carlo (VMC) and the OneShiftOnly optimizer implemented in QMCPACK. The Jastrow factors were optimized for the adsorbed systems at Γ-point only and the element-specific parameters from these optimizations were used in the reference systems (M@Gr and H 2 ) for consistency and efficiency. Note that we also invoked the DLA scheme which removes the bias of the Jastrow factor on the energy in the limit of small time steps. 26 We used time-steps of 0.03 a.u. for all systems in the fixed-phase DMC calculations and we tested additional time-steps of 0.06, 0.01, 0.005 a.u. See SM for a full account of DMC time-step convergence tests. In addition, we used 21,000 total walkers across 150 GPU accelerators. Workflows were partially automated using NEXUS. 35

III. RESULTS
In the following we report the most accurate hydrogen adsorption energies available to date for metal decorated graphene sheets using fixed-phase DMC and establish the existence of a Kubas-bound chemisorption minimum on Ca@Gr. We also report the adsorption energy of hydrogen on pristine graphene with half the stochastic uncertainty from previous work. 5 In Section III A, we report the DMC adsorption energies for fixed structures found from the screening of hydrogen molecules on M@Gr systems. 2 In Section III B we uncover the difference between chemisorption and physisorption of H 2 on Ca@Gr with DMC in order to understand the reliability of DFT methods for Kubas-type binding. In Section III C we provide further heuristic insights on the effect of electron delocalization on the adsorption energy of hydrogen molecules in Kubas-type bonding using the Hubbard U method in DFT.
A. Reference H 2 adsorption energies on Li, Na, Ca, and K decorated graphene The DMC and PBE+D3 predictions of H 2 adsorption on pristine graphene and adatom decorated graphene are reported in Table I and shown in Fig. 2. In addition, we report a more precise prediction of the hydrogen adsorption energy on pristine graphene (−20 ± 3 meV), which is consistent with our previous result (−24 ± 11 meV). 5 The DMC reference adsorption energies allow us to assess the PBE+D3 adsorption energies, keeping in mind that PBE+D3 geometries are used throughout. PBE+D3 overbinds the systems we report by 9 meV in 3H 2 +Li@Gr, 8 meV in 3H 2 +Na@Gr, 17 meV in 4H 2 +K@Gr, and 74 meV in 4H 2 +Ca@Gr, per H 2 molecule with respect to the DMC refer-  ences (outside of the 1σ stochastic error bars in DMC). Therefore, the performance of PBE+D3 across these relatively similar materials is variable, indicating that PBE+D3 and indeed similar DFT methods may not accurately rank different materials for hydrogen storage on a large-scale. The large discrepancy of 70 meV between PBE+D3 and DMC for 4H 2 +Ca@Gr is particularly noteworthy. In our previous work at the DFT level, Ca@Gr was found to be the most promising material among group 1 and group 2 adatom decorated graphene sheets for hydrogen storage. This was partly due to the favourable adsorption energy predicted using PBE+D3 and also thanks to the tuneable mechanism that underpins this binding. More specifically, it has been shown that a weak Kubas-type covalent bonding can exist between hydrogen and a metal atom with partial d-state occupation. We showed using DFT that this form of binding can be tuned with experimentally accessible controls such as external electric fields and substrates supporting graphene. Here, we find that PBE+D3 significantly overbinds 4H 2 +Ca@Gr and according to DMC, Ca@Gr and K@Gr are the weakest adsorbers of hydrogen among the four materials computed in this work. In order to understand whether Kubas-type binding of hydrogen is feasible at all on Ca@Gr we have to consider the relative energy difference between the Kubas-bound 4H 2 +Ca@Gr system and a corresponding physisorption structure on Ca@Gr. This is addressed in the following section.

B. Is H 2 Kubas-bound or physisorbed on Ca@Gr?
We have established that the PBE+D3 absolute adsorption energy for the (chemisorbed) Kubas-bound 4H 2 +Ca@Gr is overestimated by 74 meV with respect to DMC. Furthermore, the −113 ± 2 meV hydrogen adsorption energy predicted by DMC is far outside the favourable window of −200 to −400 meV for hydrogen storage. 2,3 However, it is nonetheless important to establish if the Kubas-type mechanism of binding predicted at the DFT level for Ca@Gr is feasible since we previously showed that such a covalent form of interaction has potential for being tuned towards more favourable hydrogen adsorption energies. To this end, our goal is to compare the thermodynamic binding preference between the Kubas configuration and a corresponding physisorption complex. We performed PBE+D3 geometry optimizations for several different starting positions of 4H 2 molecules centered on Ca@Gr. The initial structures included H 2 in upright and flat orientations relative to the graphene sheet and Ca-H 2 distances of ∼ 3.5 Å. The most favourable physisorption structure is shown in Fig. 3 and its adsorption energy is −81 meV per H 2 molecule with PBE+D3. Using the same physisorption configuration in DMC and a time-step of 0.03 a.u., the physisorption energy is −31 ± 3 meV. Therefore, PBE+D3 overestimates both the chemisorption and physisorption binding energies, but their relative order is preserved in DMC indicating that Kubas-type binding is thermodynamically feasible. DMC predicts that the chemisorption state is favoured over physisorption by ca. 80 meV at 0 K. As such, there is renewed potential for Kubastype interactions to be exploited as a tuneable binding mechanism for storage of H 2 in metal decorated materials.
The considerable difference between the Kubas chemisorption and physisorption energy of 4H 2 +Ca@Gr can be understood by noting the different electronic states that are the highest occupied molecular orbitals (HOMOs) in each configuration, as shown in Fig. 4. Given that the Slater determinant of the DMC wavefunction is initialized by LDA orbitals (and the nodal surface remains fixed), it is evident from Fig. 4 that the chemisorption and physisorption states are starkly different in electronic configuration as well as geometry. Therefore, it is important to consider how the electronic configuration of the chemisorption state is affected by different electronic structure methods. To this end, we gauge the effect of electron localization on the electronic structure and adsorption energy of 4H 2 +Ca@Gr in Section III C.

C. A heuristic test of electron localization using Hubbard U
Previous DMC and coupled cluster (with single, double, and perturbative triple excitations) predictions of a Ca + -4H 2 isolated gas cluster showed that Kubas-type binding is not thermodynamically favoured relative to physisorption and that DFT methods strongly favour Kubas-type binding incorrectly. 13 In this work, we also computed the Ca + -4H 2 gas cluster using our DMC protocols as well as using LDA and hybrid B3LYP. Our results are consistent with previous works and we report the findings in Appendix Section B for the interested reader. This qualitatively wrong finding from DFT can be attributed to the delocalization error in exchangecorrelation functionals, such that the stabilization from the orbital overlap of the 3d-state of Ca + and the H 2 1σ * orbitals is overestimated. Importantly, this was found to be the case even using a hybrid density approximation such as B3LYP, which partially corrects for the delocalization error by including a fraction of exact exchange.  In the materials we consider, the metal adatoms are at least partially oxidized by the graphene sheet, making the whole system metallic 2 and thus setting them apart from the gascluster system computed with reference methods in the past. To understand the impact of increasing the electron localization on the chemisorption and physisorption minima of 4H 2 +Ca@Gr, we performed heuristic PBE+D3 calculations with Hubbard U corrections for U values between 1 to 6 eV. In doing so, we find that the physisorption energy is relatively unperturbed for different U values, ranging from −82 to −86 meV, as can be seen from Fig. 3 and Table II. Whereas, the chemisorption minimum is strongly affected by electron localization, weakening with increasing U value, by up to −114 meV at U=6 eV (−71 meV after geometry relaxation). This can be understood by observing that the HOMO for the Kubas-bound system is dominated by the 3d state of Ca overlapping with the anti-bonding 1σ * orbitals of H 2 molecules, while the physisorption HOMO is dominated by the Ca 4s state, as can be seen in Fig. 4. Importantly, the HOMO of the chemisorption state is dominated by the 3d state of Ca even up to U=6 eV as can be seen in the projected density of states (PDOS) of the chemisorption system at different U values in Fig. 5. It can also be seen that the Ca 3d state is less occupied for U=6 eV than for U=1 eV. As such, the application of the U value to the d-state of Ca affects more strongly the binding energy of the Kubas-bound chemisorption structure. While this explains the sensitivity of the chemisorption energy to electron localization, it also suggests that different exchange-correlation functionals yield the same order of occupied electronic states and similar nodal surface for a given geometry. As fixed-phase DMC is constrained only according to the nodal surface, the implication is therefore that the DMC chemisorption energy is unaffected by different input DFT orbitals -as is also shown in the Ca + +4H 2 gas cluster.
Separately, we also fully relaxed the geometries at each U value and, as can be seen from the adsorption energies in parenthesis in Table II, we found that the structures do not change noticeably with the exception of PBE+D3+U where U is 5 and 6 eV. Starting from the Kubas-bound structure with high values of U results in re-configuration to a physisorption state and the bond length of H 2 molecules shorten. The H-H bond length at U of 5 and 6 eV is 0.76 Å, and the distance to the Ca adatom is ∼ 2.56 Å and ∼ 3.04 Å, respectively. In all other cases, the minimal change in structure and adsorption energy upon relaxation with different U values also suggests the use of fixed PBE+D3 geometries in DMC does not bias the results. Overall, this heuristic demonstration on the impact of localizing the electron density indicates that the Kubas-type interaction energy is sensitive to the DFT approximation used and thus suggests that such interactions are challenging to accurately predict at the DFT level.

IV. CONCLUSION
We predicted the most accurate molecular hydrogen adsorption energies available to date on pristine graphene and metal adatom decorated graphene sheets using DMC. The DMC adsorption energy of H 2 on pristine graphene is −20±3 meV. This result is consistent with what was predicted with DMC in the past, 5 but thanks to algorithmic developments and better computational efficiency, we achieved half the stochastic error here. Going beyond pristine graphene, we show that among Li, Na, K, and Ca adatoms on graphene, Li facilitates the strongest binding of molecular hydrogen, with an adsorption energy of −172 ± 6 meV per H 2 molecule. Furthermore, reference DMC predicts that Ca@Gr is the weakest adsorber of H 2 in this subset of materials -in stark contrast to previous DFT predictions. The broader implication of this is that a widely-used DFT method such as PBE+D3 cannot accurately rank hydrogen adsorption within a small subset of materials with different adatoms and therefore cannot be expected to deliver very accurate predictions in large-scale materials screening. However, we additionally computed a physisorption minimum for 4H 2 on Ca@Gr with PBE+D3 and DMC and ascertained that a chemisorption minimum is thermodynamically favoured. We have previously shown that the chemisorption minimum is underpinned by Kubas-type covalent binding and is therefore modifiable by external controls such as electric fields and substrate materials supporting the graphene sheet. As a result, the confirmation from DMC that this form of binding is thermodynamically feasible provides support for further work on exploiting Kubas-type interactions to boost the hydrogen adsorption energy for hydrogen storage applications.
Appendix A: Sampling the Brillouin zone in DMC using information from DFT The k-point convergence of the adsorption energy on Li, Na, and K decorated graphene is achieved with a single kpoint with respect to a fully converged Γ centred 5 × 5 × 1 k-point grid, as shown in Table III. The single k-point calculations were performed non-self consistently using the fully converged charge densities with LDA and VASP v5.4.4. Although the adsorption energy cannot be obtained directly from a non-self consistent calculation using Quantum Espresso, we found that the adsorption energies are in agreement between the two codes at these k-point grids and that the same convergence is achieved. Using this information, we chose the K point to produce orbitals for the DMC adsorption energy computations of 3H 2 +Li@Gr, 3H 2 +Na@Gr, and 4H 2 +K@Gr. TABLE III. Adsorption energies per H 2 molecule (in meV) from LDA using fully converged 3 × 3 × 1 and 5 × 5 × 1 k-point grids and non-self consistently (using the converged charge density) at the Γ and K points.
The 4H 2 +Ca@Gr system has a more complex convergence with respect to k-point sampling and we found that the adsorption energy at each (non-self consistently computed) k-point varies significantly from the fully converged adsorption energy. This necessitates the twist-averaging of DMC energies with orbitals obtained at different k-points. In order to find the most accurate (and feasible) grid for twist-averaging we analysed the results of several k-point grids and the electron occupancies across a dense 15×15×1 grid for the three configurations: Kubas-bound 4H 2 +Ca@Gr, physisorbed 4H 2 +Ca@Gr, and the Ca@Gr substrate. The electron occupancy at different points in reciprocal space in the first Brillouin zone can be seen in Fig. 6. In DFT, non-integer electron occupations are possible thanks to smearing across a set of k-points, however, in our real-space QMC simulations only integer electron occupations are allowed. As a result, the most suitable kpoint grids at the DFT level correspond to those that avoid non-integer occupations and achieve good convergence for the adsorption energy. The 3 × 3 × 1 grid that is shifted by (1/2, 1/2, 0) from the Γ point centre is well-converged for the adsorption energy and has mostly integer occupations at each k-point, as can be seen from Table IV and Fig. 6. Out of the 9 k-points, there are only two k-points with 1/2 electron each at the LDA level, and this is the case only for the physisorbed 4H 2 +Ca@Gr and Ca@Gr substrate systems. The two k-points with 1/2 electron occupation are equivalent in energy and as such, we can assign one full electron to only one of the k-points and 1/2 an electron less on the other k-point. In doing so, the correct total number of electrons and total energy is maintained when averaged. We report the occupancies from LDA and those used at the QMC level in Table IV.   TABLE IV. The adsorption energy and electron population from the 3 × 3 × 1 off-centred grid and at the corresponding separate k-points. The adsorption energy at each individual k-point is computed using NSCF calculations based on the converged charge density and the corresponding integer number of electrons in each row. The adsorption energy is reported for the Kubas-bound chemisorbed 4H 2 +Ca@Gr system. Ca is also found to be partially occupied with some methods, enabling a Kubas-type binding. 13 We have verified this using ccECP pseudopotentials and the QMCPACK code in this work, see Fig. 7. The gas cluster geometry is obtained from a constrained B3LYP geometry optimization, using ORCA and def2 [3/4] basis set extrapolations. The constraint was to keep the square planar orientation of the complex with Ca at the centre, while the H-H bond length is allowed to relax. In the work of Bajdich et al., the H-H bond length was kept fixed for different Ca-4H 2 separation distances instead. However, Purwanto et al. verified that the impact of flexible H-H bonds is small while predicting the gas cluster with auxilliary field QMC. 14 We used Quantum Espresso to produce LDA and B3LYP orbitals for DMC with 600 Ry planewave cut-off and the gas cluster in a unit cell box (15 × 15 × 15) Å. We find that the B3LYP and LDA orbitals yield near-indistinguishable interaction energies across the separation distances where physisorption and Kubas-type binding may occur. In addition, the DMC total energy with LDA and B3LYP orbitals is statistically indistinguishable, e.g. at a separation distance of 2.2 Å the DMC total energy with LDA orbitals is −1119.412 ± 4 eV while it is −1119.408 ± 3 with B3LYP. Moreover, we confirm that DMC does not clearly favour Kubas-type binding in the gas cluster and this is independent of the exchange-correlation functional used to obtain orbitals. The main effect of the underlying orbitals (B3LYP and LDA) is to shift the distance at which the HOMO switches from the 3d orbital of Ca (corresponding to the first minimum) to the 4s orbital (corresponding to the second minimum), as can be seen from Fig. 7. More specifically, the DMC-LDA HOMO at 2.8 Å is mostly Ca-3d while DMC-B3LYP HOMO at this distance has Ca-4s character. However, near the interaction energy minima, LDA and B3LYP orbitals are consistent and therefore the DMC energies are in good agreement. Our calculations are therefore consistent with previous QMC works. 13

S1. TIME-STEP CONVERGENCE IN DMC
The H 2 +Gr system was previously computed using the CASINO code with Trail and Needs pseudopotentials 1,2 and a time-step of 0.015 a.u. in DMC, based on LDA orbitals.
In Table S1 we report additional calculations from CASINO and QMCPACK at different time-steps, demonstrating agreement between codes and time-steps. It can be seen from Table S2 that there is only a small time-step dependence of the DMC adsorption energies using 0.03 and 0.06 a.u. time-steps. We computed the chemisorbed 4H 2 +Ca@Gr system using 0.005 a.u. time-step, resulting in an adsorption energy of −111±5 meV, which is within the 1σ stochastic error of the result from 0.03 a.u. in DFT, which are implemented for non-spin polarized systems in QMCPACK and Quantum Espresso, respectively. Our production calculations are spin-polarized at the DFT and QMC level to properly represent the electronic structures that result from metal adatoms adsorbing on graphene. As such, we estimate two-body finite effects using the KZK method in a non-spin polarized calculation of the 4H 2 +Ca@Gr adsorption energy. We find the difference between KZK and LDA to be 3 meV in the adsorption energy per H 2 molecule.
We consider this to be near-negligible and within the stochastic error that we obtain for QMC adsorption energies. Note also that the difference between spin polarized and nonspin polarized adsorption energy for this system is ca. 30 meV at the LDA level and therefore it is important to perform spin-polarized simulations.

SYSTEMS
The band structure has been computed non-self consistently for the Γ-M-K path using the charge density from a 15 × 15 × 1 k-point mesh (computed self-consistently). See