Physi-Sorption of H₂ on Pure and Boron–Doped Graphene Monolayers: A Dispersion–Corrected DFT Study

Abstract

High-surface-area carbons are of interest as potential candidates to store H₂ for fuel–cell power applications. Earlier work has been ambiguous and inconclusive on the effect of boron doping on H₂ binding energy. Here, we describe a systematic dispersion–corrected density functional theory study to evaluate the effect of boron doping. We observe some enhancement in H₂ binding, due to the presence of a defect, such as terminal hydrogen or distortion from planarity, introduced by the inclusion of boron into a graphene ring, which creates hydrogen adsorption sites with slightly increased binding energy. The increase is from −5 kJ/mol H₂ for the pure carbon matrix to −7 kJ/mol H₂ for the boron–doped system with the boron content of ~7%. The H₂ binding sites have little direct interaction with boron. However, the largest enhancement in physi-sorption energy is seen for systems, where H₂ is confined between layers at a distance of about 7 Å, where the H₂ binding nearly doubles to −11 kJ/mol H₂. These findings suggest that interplanar nanoconfinement might be more effective in enhancing H₂ binding. Smaller coronene model is shown to be beneficial for understanding the dependence of interaction energy on the structural configurations and preferential H₂ binding sites.

1. Introduction

Among hydrogen storage materials, activated carbons, graphene, carbon nanotubes, and metal-organic frameworks are of great interest due to a combination of large surface area, low mass density, and high structural stability [1]. However, the typical H₂ binding energy of these high-surface-area carbon materials is weak (–5.1 kJ/mol H₂) [2]. The weak binding energy is attributed to physi-sorption dominated by van der Waals attraction, such that adsorbed molecules of H₂ do not remain attached at ambient temperatures and require liquid nitrogen temperatures to achieve adsorption [3,4]. For practical applications, such as H₂ fuel–cell electric vehicles, the operating temperatures are much higher, ranging from −40 to 85 °C, and require binding energies on the order of 15 to 25 kJ/mol for pressures ranging from 5 to 100 bar [5,6]. Unsaturated metal centers (e.g., Mg and Ni in metal–organic frameworks) provide an effective approach for enhancing the strength of binding between H₂ and nanostructured porous materials [7,8]. Other approaches have involved templating porous carbons that can provide impressive H₂ sorption properties and gravimetric and volumetric H₂ storage capacities of approximately 5 wt% and 35 g/L, respectively, at 77 K and 20 bar [9]. Confinement effects have further been investigated in H₂ binding in potassium-graphite-intercalation compound KC₂₄, in which H₂ is found to adsorb with an enthalpy of 8.4 kJ/mol H₂ between layers of graphite with an interplanar spacing of 3.35 to 5.4 Å [10]. Similarly, graphite nanofibers have been studied as potential materials with enhanced storage capacity due to their ability to capture H₂ molecules between layers, being estimated to be ~3.4 Å apart. Ahn et al. [11] reported that the performance of graphite nanofibers is not superior to other activated carbon materials, even though other reports have given ranges of storage between ~3 wt% and up to 13 wt% [12,13,14]. Calculations of H₂ adsorption on slit pore models of graphite nanofibers by Wang et al. [15,16] also suggest that although the pores show increased adsorption relative to modeled single-wall carbon nanotubes, slit pores with a distance of 9 Å had the highest adsorption capacity. In addition, doping carbons with electron–deficient boron (B) has been investigated because of the lower mass of B compared to metal. Substituting B for a metal provides a potential approach to meeting gravimetric density targets for storing hydrogen onboard a fuel–cell electric vehicle [17]. In a pioneering computational study [18], doping fullerene C₃₆ with B was predicted to lead to a significant increase in the hydrogen adsorption energy to −19.2 kJ/mol H₂. The authors of that study suggested that the increase in binding energy is due to a partial charge transfer from the sigma-bond of an H₂ molecule to the free localized p_z–orbital of a B atom. However, a notable feature is the intrinsic non-planarity of the fullerene that induces distortion of the B center, allowing for better matching of the orbitals for enhanced interactions. In contrast, Zhou et al. [19] employed periodic one-dimensional calculations on carbon nanotubes using density functional theory (DFT) with a generalized gradient approximation (GGA) functional and supercells that range from 64 to 80 atoms and reported H₂ binding and pristine carbon material binding energies that ranged between −1.8 and −10.3 kJ/mol H₂, depending on the position of the H₂. Furthermore, they observed that B doping reduces the binding energy by between 0.7 and 3.6 kJ/mol H₂ and suggested that doping with electron–deficient B disrupts the conjugated six–membered ring electronic system. They also suggested that the polarity induced in the nanotubes by doping might play a role in the adsorption of the H₂ molecule, resulting in decreased binding energy. More recently, a theoretical study on the adsorption of various gas molecules on graphene that was doped with B [20] reported energies of adsorption of −1.35 kJ/mol H₂ using GGA, acknowledging that, for non-covalent interactions, the consequent neglect of long–range van der Waals (vdW) interactions, does not accurately describe the physi-sorption of H₂. Computational work by Sha et al. on H₂ storage in bulk BC₃ material shows significant overestimation of the binding energies and spontaneous dissociation of H₂ when local density approximation (LDA) functionals were used [21]. It has been established that LDA and GGA functionals are insufficient for capturing the vdW interactions of H₂ with doped carbon materials [19,20]. These shortcomings have led to the development of the vdW–DF correlation functionals (vdW–DF1 [22,23] and vdW–DF2 [24]), which are parameterized to capture the vdW forces without any empirical corrections to the nonlocal dispersion interactions for a variety of systems, including graphene along with other carbon–based polymeric materials. Here, the nonlocal correlations are taken into account through the density–density interaction term [22]. Klimes et al. [25] showed that the accuracy of these functionals can be dramatically improved in dispersion bonded complexes. Benchmark studies of the vdW–DF2 functional with revPBE [26] exchange on H₂ adsorption on coronene and graphene surfaces show physi-sorption energies for H₂ on the hollow site of coronene (–5.7 kJ/mol H₂) and graphene surfaces (–5.9 kJ/mol H₂) to be in good agreement with experimental values [27]. A benchmark investigation of the physi-sorption of H₂ on pure benzene and graphene shows that the revPBE based vdW–DF2 method with nonlocal dispersion is in good agreement with high–level coupled cluster (CCSD(T)) calculations and experimental values for both benzene and graphene [28].

Being encouraged by theoretical predictions of enhanced binding of H₂ to B-doped carbon materials, several experimental studies have been undertaken to dope these materials with boron. Some of these experimental approaches targeted the increased heteroatom substitution, which is the nature of the dopant in frameworks with high surface area and tunable porosity. Wang et al. [29] used the substitution reactions to prepare B-doped microporous carbons. They reported that B doping led to an increase in the surface area and the isosteric heat of adsorption [Q_st] for H₂ compared to the undoped carbon (i.e., an increase from ca. −6 kJ/mol to −10 kJ/mol at higher coverage). Chung et al. [30] developed a new class of polymeric B-containing precursors that were thermally transformed to microporous B/C materials with a boron content of 7.2% and surface area of 780 m²/g. The substitutional p–type B dopant was shown to polarize the C surface, resulting in higher Q_st of −10.8 kJ/mol for 0.62 wt % at temperatures of 77 to 87 K and a pressure of 1.2 bar for H₂ binding. Tour and coworkers [31] reported a bottom-up, solution–phase technique for the preparation of pristine and heteroatom-substituted carbon scaffolds that show high surface areas and enhanced H₂ physi-sorption capacities relative to the pure carbon scaffolds. This approach involved heating chorine-containing small organic molecules (precursors) with metallic sodium at reflux in high-boiling solvents, being subsequently followed by the addition of the heterotopic electrophiles to the mixture for dopant incorporation. The substituted carbon scaffolds enriched with 7.3% boron were observed to have higher surface area of around 900 m²/g and enhanced reversible hydrogen physi-sorption capacity (–8.6 kJ/mol at zero coverage). Jeong et al. [32] developed a method for preparing a broad range of porous B-doped carbon materials by using B precursors containing inorganic additives. They reported the chemical structure of the doped sheets changed from a disordered (less π–conjugated) state with a B-puckered configuration at 600 to 800 °C to an ordered (highly π–conjugated) state with a planar configuration at 1500 °C. The planar graphitic layers were reported to accommodate <3% of B content, while the amorphous, BC_x-like materials exhibited much higher surface areas (500 to 800 m²/g) and 12% B content with electron–deficient B moieties, resulting in enhanced H₂ binding energies (–12 to −20 kJ/mol at higher coverage). In a recent joint experimental and computational effort [33], a site–specific selective incorporation of B into the graphene sublattices was demonstrated to form a tunable band gap controlled by the dopant concentration. Shcherban et al. [34] reported a nanocasting approach to prepare B-doped CMK-3 using SBA-15 as an exo-template. They reported no increase in H₂ binding energy as a result of incorporating B and a decrease in the overall capacity of H₂ adsorption. In 2012, Bult et al. [35] reported enhanced hydrogen binding energy for a high-surface-area, B-doped carbon material made by chemical vapor deposition of a B-doped carbon (BCx) layer onto CM-Tec activated carbon. The enhancement was attributed to increased interlayer spacing of the material (~4 Å) and higher concentration of B sites.

Table 1. Binding energy of H₂ (kJ/mol H₂) on carbon and boron–doped carbon materials. The * denotes calculated values.

System	Binding Energy of H2	Reference
Pure Graphene	−5.1 kJ/mol H2	[2,3]
Microporous B/C	−10.8 kJ/mol (0.62 wt %)	[30]
B doped Carbon Scaffolds	−8.6 kJ/mol (zero coverage)	[31]
Porous BCx	−12 to −20 kJ/mol (high coverage)	[32]
B doped CMK-3	−5 kJ/mol H2	[34]
B doped C	NA	[36] (increased capacity to 5 wt% at 298 K)
C36 Fullerenes	−19.2 kJ/mol H2	[18] *
C–Nanotubes	Reduced from −10.3 kJ/mol H2 in pristine	[19] *
B doped Graphene	−1.35 kJ/mol H2	[20] *
Pure Coronene	−5.7 kJ/mol H2	[27] *
Pure Graphene	−5.9 kJ/mol H2	[27] *
Pure Graphene	−5.5 kJ/mol H2, −4.6 kJ/mol H2	[28] *

provides a summary of the previously measured and calculated binding energies.

Experimental reports of B-containing carbon materials for H₂ storage are relatively limited when compared to computational studies. The highly oxophilic nature of boron makes synthesis of these materials with precise control and purity very challenging. In addition, experimental data for H₂ capacities and binding energy show significant discrepancies, potentially due to different measurement protocols and a lack of a standard method of measurement, as well as wide range of the materials and the impurities they contain. Thus, reliable theoretical methods for predicting H₂ physi-sorption energies can be helpful in rationally designing these storage materials and achieving improved structure–activity relationships.

In this study, we report a detailed investigation of H₂ adsorption on a coronene model and a periodic graphene monolayer and bilayer, as well as their boron–doped variations. We also investigate the effect of B-doping on graphene at different substitutional sites and the effect of interlayer spacing on the binding properties of pure and doped graphene surface. We found that explicit inclusion of the vdW dispersion correction through nonlocal correlation parameterized by density and its gradient is key to capture the physi-sorption of H₂. The accuracy of these binding energies is quite sensitive to the optimal density gradient correction in the vdW kernel. Using improved DFT exchange with the dispersion corrected correlation functional provides excellent agreement with experiment [2]. The coronene model was found to adequately describe the binding of H₂ and it can be used both for exploratory and detailed electronic structure calculations of the interacting sites in lieu of the much more expensive periodic calculations. The adsorption properties of the doped coronene are also analyzed using quantum molecular descriptors, such as charge and bond order analyses, as a function of the position of B-atom substitution and the H₂ adsorption site.

2. Methods

Quantum mechanical calculations for the coronene (C₂₄H₁₂) and the doped coronene were done with the NWChem computational package [37] while using DFT with the hybrid B3LYP (a = 20) functional [38,39]. The results also were compared against the GGA functionals PBE0 [40] (a = 25) and PBE [41] (a = 0). Here, parameter a is the fraction of Hartree–Fock (HF) exchange in the Exchange–Correlation (XC) functional, given as:

$$E_{\mathrm{XC}}=a{E}_X^{HF}+\left(1-a\right)E_X^{GGA}+a{E}_C^{GGA}$$

(1)

For all calculations, we used the 6-311G++** basis set [42], and dispersion interaction was included within the framework developed by Grimme (D3) with the semiempirical correction [43]. Vibrational analysis was done to confirm the minima at the potential energy surface, at the rigid rotor-harmonic oscillator approximation, as well as to provide zero–point correction to the electronic energies. As we did not observe any significant change in the energies with the vibrational corrections, we only report electronic energies in this study for consistency with the solid state calculations. All the calculations were done in the gas phase, at 0 K temperature.

Three-dimensional periodic system calculations were done using the VASP code [44]. To model the H₂ adsorption into the graphene surface, we selected a rectangular supercell structure with the cell size of 12.3 × 12.3 Å separated by the vacuum layer with thickness of 15 Å and 7 Å between two monolayers along the z–axis. Vacuum layer with thickness of 15 Å behaves as a two-dimensional material, with the presence of the two monolayers having no impact on the energetics. The unit cell contained 50 carbon atoms. We optimized the structure using the DFT PBE [41] level of theory with projected augmented wave potential [45], with a cutoff of 875 eV. The total energy of the system was minimized with respect to the internal coordinates. The energies of self-consistent calculations were converged to the 10⁻⁶ eV/cell. A Monkhorst–Pack grid of 36 k–points (4 × 4 × 4 mesh) and 52 k–points was used (4 × 4 × 6 mesh) to sample the Brillouin zone. Both of the meshes yielded similar results. The dispersion corrections were added using the vdW–DF2 [24] nonlocal correlation functional with its semilocal PW86 exchange [41] replaced with that of GGA PBE [41], revPBE [26], optPBE [25], and RPBE [46] functionals as implemented in VASP. The performance of vdW–DF2 functional at PBE [41] exchange level was compared to the erlier vdW–DF1 [23] version and other DFT–D3 methods (with zero and Becke-Johnson damping) [43,47] with semiempirical dispersion corrections. The vdW–DF2 functional uses a gradient correction in its parameterized vdW kernel in contrast to the older vdW–DF1 version [24]. The effect of the interplanar spacing between two graphene monolayers was tested using both the vdW–DF1 [23] and vdW–DF2 [24] correlation functionals with different DFT exchange contributions [25,41,46].

3. Results

3.1. Coronene Model Calculations

Coronene is often used as a model system for a graphene surface. It is a polycyclic aromatic hydrocarbon that contains six benzene rings (Figure 1). The smaller size of coronene, allows for detailed calculations and the rapid evaluation of variety of doped structures, while maintaining the periodic aromaticity of graphene. The approach taken in this study was to use the smaller model for evaluation of the doping effects on the H₂ binding and subsequently refine the calculations for the most promising systems with periodic calculations.

The binding energy ($\Delta E_b$) was calculated as the difference between the electronic energy of the molecule with and without H₂; that is

$$\Delta E_b=E_{system+n{H}_2}-\left[E_{system}+n{E}_{H_2}\right]$$

(2)

where n is the number of hydrogen molecules adsorbed. We used the electronic energies to evaluate the binding energies of the materials, for consistency with the solid–state calculations. We do not expect significant thermal contributions since H₂ physi-sorption occurs at low temperatures.

3.1.1. Benchmarking Level of Theory for H₂ Binding to Coronene

Validation of the calculations on coronene as a model system show that adding the semiempirical Grimme dispersion correction to the DFT electronic energies significantly affects the energies for H₂ physi-sorption.

Table 2. Binding energy ($\Delta E_b$ in kJ/mol H₂) for H₂ physi-sorption on coronene using different XC functionals.

XC Functional	$\mathbf{\Delta }{\mathit{E}}_{\mathit{b}}$	$\mathbf{\Delta }{\mathit{E}}_{\mathit{b}}$ (D3 Corrected)
PBE	−0.3	−9.8
PBE0	−0.7	−7.3
B3LYP	3.1	−6.2
vdW–DF2 (solid state) [27]		−5.7

presents the binding energies for coronene with and without the DFT–D3 dispersion correction for PBE, PBE0, and B3LYP XC functionals. The lack of dispersion correction in all cases grossly underestimates the electronic binding energies, with B3LYP predicting an endothermic reaction. The addition of the Grimme D3 dispersion results in the increase in the binding of about 8 to 9 kJ/mol H₂, depending on the functional used, predicting binding energies in the range of −6.2 to −9.8 kJ/mol H₂. The binding energies calculated at B3LYP level with D3 semiempirical correction is in excellent agreement with the calculated [27] H₂ binding on coronene sheets. For all additional studies on the coronene model, we used the D3 dispersion–corrected B3LYP functional. All of the structures and absolute energies for the molecular models are provided in the Supplementary Material (Pages p8–p14).

3.1.2. Effect of Boron Substitution on H₂ Binding

We considered several plausible configurations, as indicated in Figure 1, for boron atoms occupying sites on the central ring of the coronene model. All structures were doped with either one or two B atoms. In the systems where only one B atom was introduced, the charge was preserved by terminating an adjacent carbon atom with hydrogen. The nomenclature of the systems refers to the position of the dopant on the central ring of the coronene. For the ortho position, the two dopant atoms are adjacent to each other. In the meta position, there is one C atom between the two dopants, while para means that there are two C atoms between the two dopant atoms.

Doping with one B atom enhances H₂ physi-sorption by 1.4 kJ/mol H₂ to −7.6 kJ/mol H₂, whereas the addition of a second B molecule to the central ring results in a decrease, with the strongest binding energy being only 0.4 kJ/mol H₂ higher than pristine coronene, as observed for ortho substitution. The meta and para positions slightly reduce the energy of physi-sorption in comparison to the undoped system.

Table 3. Binding energy ($\Delta E_b$ in kJ/mol H₂) of adsorbed H₂ in coronene and boron–doped models. The nomenclature introduced here will be used for shorthand notation of the systems hereon.

System	Nomenclature	$\mathbf{\Delta }{\mathit{E}}_{\mathit{b}}$
C24 H12	C	−6.2
B C23 H13	B	−7.6
B2 C22 H12	BB ortho	−6.6
	BB meta	−5.9
	BB para	−5.6

summarizes these results and the trends are discussed in detail in the following two sections. The dependence of site of physi-sorption and strength of binding energies as a function of B–doping is explained in Section 3.1.3 and Section 3.1.4, respectively.

3.1.3. Effect of Boron Doping on the Charge Distribution

We consider the distribution of charges in the central ring of the model to better understand the sites of H₂ adsorption in different models and variations in the binding energies. We report the Mulliken charges for all atoms in

Table 4. Mulliken charges for the atoms on the central ring in the coronene models. The numbering of atoms and naming of system is as shown on Figure 1 and Table 3. The substituted boron atoms are identified by a superscript ‘b’. The distorted C^d atom is identified by a superscript ‘Cd’.

Atom #	C	B	BB ortho	BB meta	BB para
1	0.10	0.01 b	−0.22	0.51 b	0.56 b
2	0.06	−0.28	0.13	−0.22	−0.21
3	0.06	0.11	0.13	0.16	−0.12
4	0.11	0.09	−0.27	−0.22	0.37 b
5	0.04	0.14	−0.03 b	0.50 b	−0.19
6	0.05	−0.54 Cd	−0.03 b	−0.56	−0.20

. There is no separation of charges in the central ring due to the symmetry of the carbon atoms in the coronene ring, resulting in H₂ interaction with the full ring, positioning in the middle of the ring, known as the ‘hollow’ site, parallel to the plane of the surface.

Single and double substitutions of boron result in either positive charge distribution on the boron centers (BB meta and BB para) or nearly neutral (B and BB ortho), and negative charge on the adjacent carbon atoms. In the singly doped B system, the puckered C atom (C_d) adjacent to boron, which has terminal H, becomes the site for H₂ physi-sorption. In the systems that were doped with two boron atoms, the H₂ molecule physi-sorbs over the BB bond (‘bridge’ site in the BB ortho) even though this site has no deficit of electron density, or in the case of the BB meta and BB para on the most neutral C atom (‘top’ site on C3) on the central ring. This suggests that H₂ prefers to interact with the coronene with sites that are neutral when no geometric distortions are present; however, the H₂ prefers to interact with the distorted site if present in spite of the excess charge that is developed at this site.

3.1.4. Effect of Boron Doping on the Bond Order

In coronene, the electrons on the carbon centers are delocalized over the central ring, thus resulting in all C–C bonds attaining partial double bond character (1.5 ≥ bond order ≥ 1.2). B–doping changes the electronic structure, thus affecting the bond order. We report the bond orders for the pure and doped coronene systems under study in

Table 5. Bond order analysis for the bonds between the atoms on the central ring in the coronene models. The numbering of atoms and naming of system is as shown in Figure 1 and Table 3. Single asterisk (*) denotes the B–C bonds, while double asterisk (**) represents the B–B bond. The bonds containing C^d atom are identified by a superscript ‘Cd’.

Bond Order	C	B	BB ortho	BB meta	BB para
1–2	1.25	1.12 *	1.24	1.01 *	0.95 *
2–3	1.21	1.30	1.06	1.16	1.43
3–4	1.25	1.14	1.25	1.17	1.00 *
4–5	1.26	1.30	1.04 *	1.01 *	1.01 *
5–6	1.23	1.06 Cd	1.24 **	1.13 *	1.42
6–1	1.26	1.00 *Cd	1.02 *	1.13 *	0.96 *

.

The analysis of the bond order in the doped coronene suggests that a singly substituted boron system maintains a bond order ~1 for the B–C bonds, which represents a single bond character. The addition of the terminating hydrogen results in all single C–C bonds for the tetra–coordinate carbon (C6). This observation is consistent with the excess electron density being accumulated in the tetrahedral carbon, which results in negative charge (see Table 4). The complete aromaticity (six partial double bonds) within the ring that was observed for the C system is partially perturbed with the formation of four single and two partial double bonds. We emphasize the charge of the B atom for this system, being less electropositive, which is attributed to the electron redistribution within the ring (see Table 4). On the other hand, the coronene system doped with two borons at ortho substitution results in an alternately single and partial double bond arrangement of bond orders. To our surprise, the B5–B6 bond in this system attains a slight partial character, suggesting that the presence of electron–rich boron centers for this system is a result of charge delocalization, being reflected as partial double bond orders. In the BB meta model, all bonds adopt nearly single bond character, consistent with accumulation of electron density on all the carbons of the central ring bound to the borons, with borons having +0.5 charge (see Table 4). However, the para substitution of the borons results in a ring with two double bonds between the C–C atoms. The H₂ binding enhancement correlates with doping that less drastically disturbs the bonding distribution (i.e., singly doped B and BB ortho), whereas a more significant redistribution of bond order and charges within the ring results in reduced binding strength (i.e., BB meta and BB para). This suggests that a complete disruption of the aromaticity in the central ring and creation of electron–deficient or electron–rich density sites leads to decreased binding energy. The binding energy is enhanced for the B-doped variations where boron is electro neutral as a result of the charge delocalization within the ring. H₂ adsorbs on the most electro neutral site in the ring when no geometrical distortions are present in the planar structure.

3.1.5. Other Possibilities for H₂ Binding

We also explored the binding of a second H₂ molecule and the possibilities of H₂ adsorption at the external rings of the coronene models, corresponding to a model where B–doping is at the edge of a graphitic monolayer. When only one B-doped site on the central ring is considered, the binding of the second H₂ is slightly weaker (by 0.3 to 0.6 kJ/mol H₂, depending on the position of the second binding event (Figure S1 in Supplementary Materials). No significant enhancement was found for H₂ binding when B-doping occurs at the boundary of the model system (Figure S2 in Supplementary Materials).

3.2. Periodic Two–Dimensional and Three–Dimensional Doped Graphene

We selected few promising models that were used to calculate the H₂ binding to a periodic two-dimensional graphene surface based on the initial screening of doped coronene model. We also considered the effects of increased boron concentration and the effects of confining the H₂ between two surfaces with an interplanar distance of 15 Å and 7 Å. Table S2 in Supplementary Materials illustrates the absolute energies for all layered structures.

3.2.1. Benchmarking Level of Theory for H₂ Binding on Graphene

The experimental binding energy for H₂ on graphene is measured at −5.1 kJ/mol H₂ [2]. The benchmarking of the level of theory used was done on pristine graphene sheets. The initial relaxation of the surface was done using the PBE XC functional, for a periodic surface with 12.3 × 12.3 × 15 Å cell size. Single point energies were calculated for this structure with different dispersion correction methods used in this work (

Table 6. Comparison of H₂ binding energy ($\Delta E_b$ in kJ/mol H₂) for different dispersion correction and density functional theory (DFT) exchange methods for interlayer spacing of 15 Å between periodic graphene sheets. The calculated values are in good agreement when the vdW–DF2 dispersion correction is used at different DFT exchange levels, with the RPBE providing the best agreement with experiment.

DFT Exchange	Dispersion Method	$\mathbf{\Delta }{\mathit{E}}_{\mathit{b}}$
PBE	no dispersion	−0.75
PBE	D3—no damping	−0.8
PBE	D3—Becke-Johnson damping	−0.75
PBE	vdW–DF1	−9.5
PBE	vdW–DF2	−6.0
revPBE	vdW–DF2	−3.8
optPBE	vdW–DF2	−4.6
RPBE	vdW–DF2	−5.2

, top four entries). The H₂ binding energy was significantly underestimated with PBE (–0.8 kJ/mol H₂) for all entries, except those using vdW–DF1 and vdW–DF2. Whereas vdW–DF1 overestimates the binding (Table 6, and Table S1 in Supplementary Materials), vdW–DF2 provides the best agreement with experiment [2] for graphene at −6.0 kJ/mol H₂. We also investigated the sensitivity of the binding energy on the DFT exchange, using vdW–DF2 with three other functionals: revPBE, optPBE, and RPBE (Table 6, bottom four entries). RPBE exchange provides even closer agreement with experiment [2] at −5.2 kJ/mol H₂ with vdW–DF2. Slightly weaker binding energy was calculated with revPBE and optPBE.

3.2.2. Effect of Boron Doping Concentration on H₂ Binding

We substituted a single carbon atom with boron for a C–to–B ratio of 49:1 (2% doping) and four carbon atoms for C–to–B ratio of 46:4 (7.4% doping) because the coronene model shows largest enhancement of binding energy for the singly doped B site in the central benzene ring (see Table 3). We find that for ~2% B–doping, the enhancement in H₂ binding energy increases to −6.6 kJ/mol H₂ from that of −5.2 kJ/mol H₂ for pure graphene (Figure 2, top and middle). From the several plausible configurations for H₂ physi-sorption sites that ranged from the ‘top’ site on puckered C (C^d) atom to the ‘bridge’ B–C and ‘hollow’ sites (Figure S3 in the Supplementary Materials), most favorable interaction was found with the C^d atom (‘top’ C^d and ‘bridge’ B–C^d sites), as also observed with the coronene model (Figure 1). At the hollow site, the binding energy is slightly lower at −5.8 kJ/mol H_2.

Previous experimental study [31] on carbon scaffolds has revealed an increase up to −8.6 kJ/mol H₂ in the physi-sorption energy with about 7.3% of boron substitution in the C network. Because the coronene model predicts largest increase when one boron is substituted per carbon ring, in the graphene monolayer with 7.4% doping we substitute one boron atom per alternate ring (Figure 2, bottom). The binding energy in this case is −7.2 kJ/mol H₂, which suggests that increasing the B concentration from 2% to 7.4% slightly increases the H₂ binding.

Because higher concentrations are not likely to be experimentally achieved with good dispersion of B sites, and due to the low enhancement in binding energy, we don’t anticipate significant benefit from exploring higher boron doping concentrations.

3.2.3. Effect of Interlayer Distance Between Graphene Monolayers on the H₂ Binding

The addition of periodicity in the direction perpendicular to the plane of the graphene sheets, allowed us to evaluate the effects of multi–layered structure on the H₂ binding. Interestingly, the position and orientation of H₂ does not change when the spacing between the surfaces is reduced from 15 Å to 7 Å, implicating long–range interactions with the second surface (Figure 3). For graphene, a bilayer separated by 7 Å, results in a significant increase of the binding energy to −10.6 kJ/mol H_2. This finding is consistent with previous calculations by Patchkovskii et al. [48], which utilize post–HF methods, giving confidence in the DFT methodology employed here.

Interestingly, a bilayer of B-doped graphene with doping concentration of 2% (C–to–B ratio of 49:1) and a separation of 7Å gave similar binding energy of (–10.3 kJ/mol H_2.) as the bilayer of undoped graphene. This observation suggests that, even though some enhancement might be gained from doping, or distortion of the surface, a much larger effect is achieved when H₂ is confined between two surfaces.

4. Discussion and Conclusion

Systematic DFT calculations for both the coronene and the two-dimensional periodic sheets of graphene systems demonstrated the importance of explicitly considering the effects of long–range van der Waals forces to the binding energies by including semiempirical (DFT–D3) and nonempirical (vdW–DF) dispersion correction, respectively. We find that using the vdW–DF2 method with the corrected density gradient term is necessary in order to accurately evaluate the H₂ physi-sorption energy for graphene. We also find that judicious selection of the underlying exchange of the vdW–DF2 functional is important, with the alternative revision of PBE functional, RPBE, showing significant improvement over the previously recommended revPBE. We benchmarked the binding energy calculations at this level of theory for graphene surface, and find that RPBE/vdW–DF2 gives the best agreement.

Our calculations show that only marginal improvements in binding energy are observed when the doping does not induce any distortion to the surface, even though earlier calculations have suggested possible enhancement to the binding energy by doping graphene with boron atoms. This is consistent with the lack of experimental evidence of increased binding that would confirm the earlier work. More significant changes to the binding energy (up to 2 kJ/mol H₂) were found in cases when the doping induces a defect, such as distortion from planarity, i.e. puckered carbon atom adjacent to boron, or hydrogenation of an unsaturated bond. This has also been observed in the H₂ interaction with BX₃:NH₃ (X = H,F,Cl) frustrated Lewis acids [49]. These types of defects may provide another avenue for exploring irregular graphene–based materials, provided they can be introduced in a controlled manner. When a defect is present, the H₂ molecule tends to adsorb at the defect site. Although the effect of surface defects can be somewhat different in the bulk material, they are still likely to affect adsorption energies and should be considered for future experiment design. In addition, the coronene calculations provide an excellent model for understanding the governing criteria for the site selection and strength for H₂ physi-sorption. Analysis of the coronene calculations show that change in the charge distribution (delocalization) and the small disturbance of the aromaticity correlate with the enhancement of binding.

The graphite monolayer models allowed for us to investigate the possible enhancement of the binding energy by increasing the boron concentration. Increasing the concentration from 2% to 7%, results in enhancement of only 0.6 kJ/mol H₂.

A notable enhancement, resulting in binding energy of −11 kJ/mol H₂ is accomplished by confinement of an H₂ molecule between two layers of graphene, with the layers being separated by a distance of 7 Å. This models a scenario in which the H₂ molecule interacts with both surfaces, resulting in a largely additive effect. Increasing the distance between the layers to 15 Å no longer benefits from this additive effect and no substantial increase in binding is observed. The addition of boron sites in the case of bilayer models shows no enhancement in the binding energies. This observation suggests that another option to increase the H₂ binding might be accomplished by understanding the effects of the confinement on the H₂ adsorption, consistent with previous findings [9,10,15,16].

Direct comparison of calculations with experiment is still difficult since it is difficult to quantify the distribution of boron in the carbon lattice in experiments [31,33], however, this work provides useful insight into both the atomistic details that may guide the preferential binding sites for H₂ in these materials, as well as the possibility of enhancing this binding by the intentional design of defect sites and confinement in the material. These approaches may be more effective in developing materials that can reach the U.S. Department of Energy targets for onboard H₂ storage.