Interactions between a H₂ Molecule and Carbon Nanostructures: A DFT Study

Abstract

On a long path of finding appropriate materials to store hydrogen, graphene and carbon nanotubes have drawn a lot of attention as potential storage materials. Their advantages lie at hand since those materials provide a large surface area (which can be used for physisorption), are cheap compared to metal hydrides, are abundant nearly everywhere, and most importantly, can increase safety to existing storage solutions. Therefore, a great variety of theoretical studies were employed to study those materials. After a benchmark study of different van-der-Waals corrections to Generalized Gradient Approximation (GGA), the present Density Functional Theory (DFT) study employs Tkatchenko–Scheffler (TS) correction to study the influence of vacancy and Stone–Wales defects in graphene on the physisorption of the hydrogen molecule. Furthermore, we investigate a large-angle (1,0) grain boundary as well as the adsorption behaviour of Penta-Octa-Penta (POP)-graphene.

1. Introduction

On the quest for environmentally-friendly energy generation, chemical energy released during formation of a water molecule out of hydrogen atoms reacting with atmospheric oxygen has been proposed as a possible route. For any practical usage, however, it is critical to solve a question of storing hydrogen. The currently available technology using high-pressure tanks is too risky for, e.g., automobile applications, and therefore alternative solutions are needed. Here, carbon-based nano-structured materials have been suggested as promising candidates for hydrogen storage [1]. In fact, it was the discoverers [2] of graphene themselves, who first showed that the novel 2D material could store hydrogen easily at cryogenic temperatures and release it again at higher temperatures [3]. Consequently, the adsorption properties of pristine graphene [4,5,6,7,8,9] were extensively studied previously both experimentally as well as theoretically. In order to tune the material response further, chemical as well as structural modifications to graphene can be introduced. These include, e.g., metal decoration [10,11,12,13] or application of perforated graphene [14,15]. The present study thus focuses on the latter area which has not been thoroughly explored yet. In doing so, we consider impact of structural defects on the local adsorption properties and try to extract simple rule-of-thumb guidelines which can be used to steer further optimization of the carbon-based materials for H${}_2$ storage.

Namely, we investigated a single vacancy and a Stone–Wales defect [16,17] as examples for point defects. As a more complex structural line defect, we considered also an energetically favorable [18] symmetric large angle grain boundary composed of (1,0) dislocations. Such defects in graphene were studied before with respect to their impact on mechanical or electronic properties [18,19,20,21,22], here we newly present their impact on the H${}_2$ adsorption capacity of graphene. Based on the results obtained from the (1,0) grain boundary we furthermore investigated Penta-Octa-Penta (POP-)graphene as a graphene modification containing no hexagons, which was theoretically predicted by Wang et al. and proposed to be a promising candidate for lithium-ion batteries [23].

All our calculations include an impact of long-range dispersion (van-der-Waals) forces as an add-on to standard density functional theory calculations using Tkatchenko–Scheffler (TS) method [24]. As this method is not (yet) widely employed, we started with a benchmark study of the hydrogen adsorption on pristine graphene, whereby we could compare our data with literature data [4], which provides a useful reference for any future studies using TS method.

Finally, in order to explain the different adsorption behavior of the studied structural defects, we implemented the TS method using the SciPy [25,26] stack, by which we can demonstrate that the changes in electronic structure result in only minor contributions to the van-der-Waals-like interactions.

2. Computational Methods

All quantum-mechanical calculations were carried out with the Vienna Ab-initio Simulation Package (VASP) [27,28]. The exchange and correlation functional was treated at generalized gradient approximation level parametrised by Perdew–Burke–Ernzerhof (PBE-GGA) [29,30]. The projector augmented wave (PAW) method [31,32] was used to describe the electron-ion interactions. The used cutoff energy of the plane-wave basis set was tested for convergence of the total energy to ≈0.1 meV/at. The tests were carried out by employing $\Gamma$-centred k-mesh sampling of the Brillouin zone for hexagonal cells (graphene), and Monkhorst-Pack [33] mesh for rectangular cells. The meshes used in the present study are summarized in

Table 1. k-meshes used for sampling the Brillouin zone.

System	k-Mesh Type	k-Mesh	Supercell
graphene	$\Gamma$-centered	$5\times 5\times 1$	$5\times 5\times 1$
single vacancy	$\Gamma$-centered	$5\times 5\times 1$	$5\times 5\times 1$
Stone–Wales defect	$\Gamma$-centered	$5\times 5\times 1$	$5\times 5\times 1$
$(1,0)$ grain boundary	Monkhorst-Pack	$7\times 10\times 1$
POP-graphene	Monkhorst-Pack	$5\times 9\times 1$

. In order to take into account also the van-der-Waals interaction, various dispersion correction methods (DFT-D2 [34], DFT-D3 [35], DFT-TS [24,36]) implemented in VASP were tested for their ability to reproduce experimental lattice parameters of graphite and graphene. The results of the benchmark study are shown in

Table A1. Benchmark study of different van-der-Waals corrections implemented in VASP in order to reproduce the experimental lattice parameters of graphite. (D2 = DFT-D2, D3 = DFT-D3, D3-BJ = DFT-D3 with Becke-Jonson (rational) damping, TS = Tkatchenko-Scheffler method, TS+IHiPart = Tkatchenko-Scheffler method with iterative Hirshfeld partitioning, MBD=Many-body dispersion energy method (MBD@rSC), dDsC = density-dependent-dispersion-correction method). IVDW refers to the setting of the van-der-Waals-method correction tag in the VASP control input file (INCAR). The last two rows show the deviation norm as well as the percentage norm. A $\Gamma$-centered $25\times 25\times 10$k-mesh was used in combination very fine FFT-grid was used (increased setting for NGXF, NGYF, NGZF parameters, and a cut-off energy of 1100 eV) in order to get accurate energies and forces.

Method	Unit	None	D2	D3	D3-BJ	TS	TS+IHiPart	MBD	dDsC
IVDW		0	10	11	12	20	21	202	4
	$c_{\exp }=6.708\mathring{\mathrm{A}}$
$c_{\mathrm{calc}}$	$\mathring{\mathrm{A}}$	7.978	6.447	6.952	6.709	6.708	6.710	6.743	6.707
$|c_{\exp }-c_{\mathrm{calc}}|$	$\mathring{\mathrm{A}}$	1.270	0.261	0.244	0.001	0.000	0.002	0.035	0.001
$|c_{\exp }-c_{\mathrm{calc}}|/c_{\exp }$	%	18.933	3.891	3.637	0.015	0.000	0.030	0.522	0.015
	$a_{\exp }=2.461\mathring{\mathrm{A}}$
$a_{\mathrm{calc}}$	$\mathring{\mathrm{A}}$	2.467	2.463	2.466	2.466	2.460	2.460	2.461	2.455
$|a_{\exp }-a_{\mathrm{calc}}|$	$\mathring{\mathrm{A}}$	0.006	0.002	0.005	0.005	0.001	0.001	0.000	0.006
$|a_{\exp }-a_{\mathrm{calc}}|/a_{\exp }$	%	0.244	0.081	0.203	0.203	0.041	0.041	0.000	0.244
$\sqrt{|a_{\exp }-a_{\mathrm{calc}}{|}^2+{|c_{\exp }-c_{\mathrm{calc}}|}^2}$	$\mathring{\mathrm{A}}$	1.270	0.261	0.244	0.005	0.001	0.002	0.035	0.006
$\sqrt{{\left|{\displaystyle \frac{a_{\exp }-a_{\mathrm{calc}}}{a_{\exp }}}\right|}^2+{\left|{\displaystyle \frac{c_{\exp }-c_{\mathrm{calc}}}{c_{\exp }}}\right|}^2}$	%	18.934	3.892	3.643	0.204	0.041	0.050	0.522	0.244

. Based on these results, all subsequent geometric optimizations were performed using the Tkatchenko-Scheffler (DFT-TS) method, performance of which was assessed already in various previous studies for carbon structures [37,38,39]. Geometric optimizations were carried out by employing a conjugate gradient algorithm with an energy convergence criterion of ${10}^{-4}$ eV (per simulation box). The energy convergence criterion for self-consistency of the electronic loop was set ${10}^{-6}$ eV (per simulation box). To minimize periodic-boundary image interactions we used a 25 Å vacuum layer along the c axis for all cells, as well as a $5\times 5\times 1$ supercell of a standard 2-atom graphene unit cell.

2.1. Graphene Adsorption Curves

In order to study the adsorption of the H${}_2$ molecule, the hydrogen molecule was placed in different alignments on three different adsorption sites (Figure 1). Directly above a carbon atom, hereafter referred to as a top site, above a bond between two carbon atoms (a bridge site) and above the centre of a hexagon (a hexagon site). In the x(y)-alignment, the H${}_2$ molecule longitudinal axis is parallel(perpendicular) to the a lattice vector, while in the v-alignment the molecule is parallel to the c lattice vector of the simulation cell. Here we note that all other papers just examine only 6 different configurations [4,40,41,42], corresponding to the three adsorption sites and x and v-alignments. Indeed, also our calculations show that the x and y alignments yield identical results, and therefore hereafter we present them as h-alignment. We also note that our setup could, in principle, provide also information as a function of rotational degree of freedom similar to, e.g., Ref. [43]. However, in the case of the adsorption maps being investigated here, such procedure becomes computationally prohibitive and hence we stick to one special alignment. For each of the resulting 9 H${}_2$ configurations, the following procedure was applied:

• Fully relax the position and geometry of the H${}_2$ molecule, while keeping the C atoms in the graphene plane fixed.
• Vary the height of the molecule with fixed geometry above the graphene sheet ranging from 2 Å to 12 Å.
• Calculate the interaction energy $E_{int,{\mathrm{graphene}+\mathrm{H}}_2}$ as:

  $$\begin{array}{c}\hfill E_{int,{\mathrm{graphene}+\mathrm{H}}_2}=E_{tot,{\mathrm{graphene}+\mathrm{H}}_2}-E_{tot,\mathrm{graphene}}-E_{tot,{\mathrm{H}}_2}\left(d_{\mathrm{case}}\right)\end{array}$$

  (1)

where $d_{\mathrm{case}}$ refers to the bond-length obtained from the full relaxation of the molecule. The energy for the hydrogen molecule $E_{tot,{\mathrm{H}}_2}\left(d_{\mathrm{case}}\right)$ with bond-length of $d_{\mathrm{case}}=0.75\AA$, which corresponds to the theoretical equilibrium bond-length, was obtained from a separate calculation. Therefore, the two H atoms of an isolated H${}_2$ molecule were separated, while for each separation the total energy $E_{tot,{\mathrm{H}}_2}\left(d_{\mathrm{case}}\right)$ was calculated. To obtain consistent results we employed DFT-TS correction, with a $\Gamma$-point only k-mesh to avoid interactions originating from periodic boundary conditions.

In the following paragraphs, $E_{ads}$ refers to the minimum of the interaction energy curve $E_{int}\left(h_{ads}\right)$, where the height above the graphene plane, $h_{ads}$, denotes the vertical distance between graphene and the center of mass of the H${}_2$ molecule.

2.2. Adsorption Maps

All adsorption maps shown in this work have been calculated with z-direction (perpendicular to the graphene plane) H${}_2$ molecule alignment, in order provide data which is comparable to literature. In order to minimize computational effort, only an irreducible wedge of the supercell was sampled for each system and subsequently unfolded to span the whole supercell (see

Table 2. Parameters of the sampling used for constructing the adsorption maps. Sampling lattice refers to the symmetry of the sampling lattice, irreducible wedge defines the directions enclosing the sampling area. Symmetry operations refer to how the irreducible wedge is unfolded in order to sample the whole supercell lattice.

System	Supercell Lattice	Sampling Lattice	Irreducible Wedge	Samples	Symmetry Operations
graphene	$a=b=12.35$ Å hexagonal	hexagonal	$\left[100\right]$, $\left[010\right]$	121
single vacancy	$a=b=12.35$ Å hexagonal	triangle	$\left[110\right]$, $\left[1\overline{2}0\right]$	66	$\underset{\left[001\right]}{\underbrace{3}}\underset{\left\{1\overline{2}0\right\}}{\underbrace{m}}$
Stone–Wales defect	$a=b=12.35$ Å hexagonal	hexagonal (rotated by 30${}^{\circ }$ with respect to the supercell lattice)	$\left[110\right]$, $\left[\overline{1}10\right]$	110	$\underset{\left\{110\right\}}{\underbrace{m}}\underset{\left\{\overline{1}10\right\}}{\underbrace{m}}$
$(1,0)$ grain boundary	$a=22.66$ Å, $b=6.54$ Å orthorhombic	rectangular	$\left[100\right]$, $\left[010\right]$	170	$\underset{\left\{100\right\}}{\underbrace{m}}\underset{\left\{010\right\}}{\underbrace{m}}$
POP-graphene	$a=11.05$ Å, $b=9.09$ Å orthorhombic	rectangular	$\left[100\right]$, $\left[010\right]$	285	$\underset{\left\{100\right\}}{\underbrace{m}}\underset{\left\{010\right\}}{\underbrace{m}}$

). For each sampling point, adsorption energy was calculated for at least 7 different molecule heights around the expected minimum, and fitted with a cubic fit in order to obtain the corresponding adsorption energy and height. Figure 3 shows clips of the calculated adsorption maps.

3. Results

3.1. Hydrogen Adsorption on Pristine Graphene

As mentioned above, we start with testing the DFT-TS method against literature data, which was obtained by employing Grimme’s D3 method (DFT-D3). Similarly to the previous works [4,8,40,41,42,44], we consider 3 high-symmetry adsorption sites (hexagon–centre of a hexagon, top–above a C atom, bridge–above a C–C bond), with both horizontal and vertical molecule alignment. The results for H${}_2$ physisorption on pristine graphene are summarized in

Table 3. Adsorption energies of the H${}_2$ molecule on the three adsorption sites in all here considered configurations (see Figure 1, x- and y-alignments yield the same behaviour and are labelled h). DFT-PBE-D3 values are taken from Ref. [4].

Configuration		DFT-PBE-TS		DFT-PBE-D3 [4]
Site	Alignment	${\mathit{E}}_{\mathbf{ads}}\left[\mathbf{meV}\right]$	${\mathit{h}}_{\mathbf{ads}}\left[\mathit{\AA }\right]$	${\mathit{E}}_{\mathbf{ads}}\left[\mathbf{meV}\right]$	${\mathit{h}}_{\mathbf{ads}}\left[\mathit{\AA }\right]$
hexagon	v	−69	3.15	−69	3.10
bridge	v	−63	3.19	−63	3.18
top	v	−62	3.19	−63	3.18
hexagon	h	−75	2.94	−67	3.10
bridge	h	−63	3.09	−58	3.16
top	h	−62	3.10	−57	3.16

and the complete adsorption curves are shown Figure 2. Our $E_{ads}$ values are in good agreement with previous works [4,42], however our work suggests the hexagon-h alignment to be the most favorable in contrast to the hexagon-v geometries reported previously. We note that also the equilibrium adsorption heights are in reasonable agreement between the DFT-TS (used in the present study) and DFT-D3 employed previously [41,42], except for the hexagon-h arrangement. Here, the DFT-TS yields somewhat tighter binding ($2.94\AA$ compared to $3.10\AA$ for the DFT-D3 method), as well as a stronger adsorption tendency. The same study was also carried out for POP-graphene, with in total nine adsorption sites (Figure A1). The corresponding results are shown in Figure A2.

3.2. Adsorption Maps and Pristine Graphene

The main focus of the present work is to investigate the role of defects on the adsorption properties. The adsorption maps were obtained for the H${}_2$ molecule in hexagon-v configuration (cf. Section 2.2). Figure 3 shows only $10\times 10$ Å${}^2$ clips of the calculated supercells, for the full adsorption maps we refer the reader to the supplementary material. The color-code and the contour lines refer to the adsorption energy and are consistent throughout the whole manuscript to allow for an easy comparison between different adsorption scenarios in various figures.

As a measure to quantify the adsorption capacity improvement compared with pristine graphene, we used the fraction of the adsorption area which exhibits more negative adsorption energy (i.e., stronger physisorption) than the minimum of graphene:

$$f=\frac{1}{A}{\int }_{A}p\left(\mathbf{r}\right)d^2\mathbf{r},p\left(\mathbf{r}\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{E}_{ads}\left(\mathbf{r}\right)\le min\left(E_{ads,\mathrm{graphene}}\right)=-68.7\mathrm{meV},\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(2)

where A is the adsorption area per supercell. This quantification is motivated by one of the main graphene advantages, namely its vast surface area, which is also relevant for other carbon nanostructures. The adsorption area is expected to correlate with the hydrogen storage capacity. Nevertheless, the area with strong adsorption capability (at $0$K, see last column in

Table 4. Summary of the adsorption characteristics: minimum and maximum values of adsorption maps for different scenarios (Figure 3, fraction of the surface (in the used supercell) with adsorption energy in the range of adsorption minimum of pristine graphene $E_{ads,\mathrm{graphene}}^{\min }+1$ meV and $E_{ads,\mathrm{graphene}}^{\min }+5$ meV, and fraction of the supercell area with adsorption energy lower than that of pristine graphene. The values of +1 meV and +5 meV, were chosen arbitrarily.

System	${\mathit{E}}_{\mathbf{ads}}\left[\mathbf{meV}\right]$		Surface [%]		Increased Surface [%]
	min	max	${\mathit{E}}_{\mathbf{ads},\mathbf{graphene}}^{\mathbf{min}}+5\mathbf{meV}$	${\mathit{E}}_{\mathbf{ads},\mathbf{graphene}}^{\mathbf{min}}+1\mathbf{meV}$	${\mathit{E}}_{\mathbf{ads},\mathbf{graphene}}^{\mathbf{min}}$
graphene	−68.7	−62.7	59.41	4.95
single vacancy	−68.9	−59.6	32.02	1.47	0.14
Stone–Wales defect	−70.2	−59.1	46.88	4.93	1.43
$(1,0)$ grain boundary	−72.6	−62.4	94.09	38.48	28.23
POP-graphene	−77.4	−64.2	100.00	71.39	49.91

) is only one of the parameters influencing the adsorption capacity. For example, external parameters such as temperature or pressure, or interaction with other H${}_2$ molecules are also influencing the overall H${}_2$ uptake, and thus we refrain from its direct relation to the hydrogen storage capacity.

As a first step, we evaluated the adsorption energy, $E_{ads}$, for the pristine graphene supercell in order to get a reference value for comparison. According to Table 3 and Figure 2 the adsorption map (Figure 3a) clearly shows that the minimum of the adsorption energy is in the center of hexagon, while the absolute value around bonds (bridge sites) and carbons atoms (top sites) exhibits a weaker interaction (difference ≈$6\mathrm{meV}/$H${}_2$).

3.3. Point Defects

Although single vacancies are hard to observe experimentally [45], they are easily tractable with DFT. It turns out that the investigated single vacancy lowers the adsorption strength in its direct adjacency with respect to the pristine-like graphene, as highlighted by the yellowish regions around the defect in Figure 3b. However, the region of influence is spatially very confined, since the adsorption energy values around second nearest neighboring hexagons already approach those of pristine graphene. Furthermore, Figure 3b suggests that the H${}_2$ molecule is less likely to be adsorbed at those hexagons which are directly adjacent to the single vacancy defect, from an energetic point of view. Furthermore, there is no improvement of the minimum adsorption energy (Table 4), the lowest value occurs along the bonds to the hexagons adjacent to the defect.

Another very common point defect in hexagonal two-dimensional materials, and thus also observed experimentally [45,46], is a rotated bond which yields the Stone–Wales defect (Figure 3c). The highest value of the adsorption energy (i.e., the weakest physisorption) of the whole cell is predicted along the heptagon bonds, although the quantitative differences compared to the highest value in pristine graphene (see Table 4) are small ($\Delta E_{ads}^{\max }\approx 2.5\mathrm{meV}$). Furthermore, we observe a narrow spatial confinement of the region with changed adsorption properties. Similarly to the vacancy, no qualitative improvement around the Stones-Wales defect can be observed as, e.g., there is no nearby region with improved adsorption behavior ($E_{ads}\le min\left(E_{ads,\mathrm{graphene}}\right)$). The only location with potentially improved adsorption are the hexagons adjacent to the pentagons, however, the quantitative improvement compared to pristine graphene (≈$1.5\mathrm{meV}$) is negligible. Interestingly the single vacancy defect weakens the minima in the adjacent hexagons, such that the surface area which lies around the graphene minimum (fourth and fifth column in Table 4) while it remains constant for the Stone–Wales defect, when comparing the values to those of pristine graphene. Nevertheless the comparison of the values of the single vacancy and Stone–Wales defect clearly shows a reduction of the adsorption capabilities.

3.4. Line Defects

Inspection of the last column of Table 4 reveals that point defects do not result in areas with lowered adsorption energy compared to pristine graphene. Moreover, the adsorption landscape is apparently more undulating as demonstrated by smaller fractions of the surface with adsorption energies close to the lowest adsorption energy. Such behaviour can even counteract physisorption, hence we continue with investigating other defects, namely the $(1,0)$ grain boundary as a representative of line defects. This is still tractable with DFT with reasonable computational demands. Due to the periodic boundary conditions, the simulation cell contains two boundaries, whose cores are $11.33\AA$ apart. Figure 3e shows a clip around the grain boundary core which itself forms a heptagon-pentagon-2 hexagon line. The minima of adsorption are located in the centers of the hexagons neighboring with the pentagons. Furthermore, the blueish color of Figure 3e, indicating less spatial confinement of the line defect, when compared to the point defects, is underpinned by the number from Table 4. For this defect already nearly one third of the area lies below the minimum graphene level.

3.5. Globally Modified Graphene

The adsorption maps of the $(1,0)$ symmetric tilt grain boundaries reveal that when measuring adsorption by the adsorption energy, a qualitative enhancement occurs around regions with pentagonal arrangement of the carbon atoms. A recently proposed globally modified graphene composed of ordered pentagons and octagons, POP-graphene, seems therefore a promising candidate for enhancing the H${}_2$ adsorption. Indeed, our calculations predict very good adsorption behavior in the pentagons area (see Figure 3e), the absolute minimum is located in the center of the octagons. With an absolute value of $-77.4\mathrm{meV}$, this minimum is the strongest of all structures investigated in the present study. Consequently, the POP-graphene shows stronger adsorption throughout the whole super cell with $\approx 50\%$ of area exhibiting lower $E_{ads}$ than $min\left(E_{ads,\mathrm{graphene}}\right)$). Moreover, the highest absolute value of $max\left(E_{\mathrm{ads},\mathrm{POP}-\mathrm{graphene}}\right)$ (weakest adsorption) is just a few meV away from the minimum of pristine graphene $\approx min\left(E_{\mathrm{ads},\mathrm{graphene}}\right)+4.5\mathrm{meV}$, resulting that the whole cell is lying below the $5\mathrm{meV}$ and over 70% of the surface below the $1\mathrm{meV}$ border respectively.

4. Discussion

In this section we try to analyze the origin of the improved adsorption of the POP-graphene as well as the apparent local adsorption enhancement cased by point defects (Figure 3). There are two possible sources: (1) change in the electron distribution within carbon-structure due to defects, that affects the polarizability; (2) pure geometry, i.e., arrangement of the C atoms. Inspecting Equation (1), the former contribution should mostly cancel out as it affects (by almost the same amount) both terms $E_{tot,\mathrm{graphene}+{\mathrm{H}}_2}$ and $E_{tot,\mathrm{graphene}}$. Therefore, the change in adsorption behavior should be linked to the changed interaction between the molecule and the graphene structure which is mainly van-der-Waals-like. This can be demonstrated, e.g., by comparing the adsorption energies the top-v configuration with $E_{ads,\mathrm{graphene}}^{\mathrm{top}-\mathrm{v}}=-62\mathrm{meV}$ (Table 3) and without $E_{ads,\mathrm{graphene}}^{\mathit{top}-\mathrm{v}}=-12\mathrm{meV}$ (pure GGA-PBE) TS correction, which agrees well with previously published data [40,42].

To bring more insight, we now analyze different contributions to the adsorption energy taking POP-graphene as an example. As the total interaction energies may be splitted up into a small PBE contribution and a dispersion part arising from the vdW corrections one can rewrite it to:

$$E_{\mathrm{int}}=E_{\mathrm{int}}^{\mathrm{PBE}}+E_{\mathrm{disp}}=\underset{\Delta E_{\mathrm{int}}^{\mathrm{PBE}}-\mathrm{Figure}4\mathrm{a}}{\underbrace{E_{\mathrm{int}}^{\mathrm{PBE}}-max\left(E_{\mathrm{int}}^{\mathrm{PBE}}\right)}}+\underset{\Delta E_{\mathrm{disp}}-\mathrm{Figure}4\mathrm{b}}{\underbrace{E_{\mathrm{disp}}-max\left(E_{\mathrm{disp}}\right)}}+C=\Delta E_{\mathrm{int}}^{\mathrm{PBE}}+\Delta E_{\mathrm{disp}}+C$$

(3)

where $C=max\left(E_{\mathrm{int}}^{\mathrm{PBE}}\right)+max\left(E_{\mathrm{disp}}\right)$ is a constant term throughout the whole simulation cell, and is found to be $C=-61.18\mathrm{meV}$ for POP-graphene. Firstly, we extracted the “adsorption maps” for pure GGA-PBE without any vdW correction. These values are shown in Figure 4a as relative values with respect to the maximum value identified on the whole adsorption surface ($\Delta E_{\mathrm{int}}^{\mathrm{PBE}}$). This way we hope to highlight the “topology” changes due to different levels of approximation. Most importantly, the adsorption surface for the pure GGA-PBE is rather flat with the difference between the maxima and minima being $min(\Delta E_{\mathrm{int}}^{\mathrm{PBE}})=-5.29\mathrm{meV}$ with the minima (darkest regions) being located in the centres of the pentagons. This is clearly different from the vdW part of the adsorption energy (Figure 4c): the adsorption minima are located in the centres of octagons and their relative depth is $min(\Delta E_{\mathrm{disp}})=-16.21\mathrm{meV}$. The adsorption surface maxima are located on the sides between octagons, same as in the case of the full adsorption surface (Figure 3e). Therefore, we ascribe the main origin of the adsorption maps topology to the vdW dispersion interaction.

As the next step we calculated dispersion energy for a hydrogen molecule for each sampling point according to Grimme’s D2 method [34].

$$\begin{array}{c}\hfill E_{disp}=-{\displaystyle \frac{1}{2}}\sum _{i=1}^N\sum _{j=1}^N\sum _{\mathbf{L}}{\displaystyle \frac{C_{6ij}}{r_{ij,\mathbf{L}}^6}}{f}_{d,6}\left(r_{ij,\mathbf{L}}^6\right)\end{array}$$

(4)

$$\begin{array}{c}\hfill f_{d,6}\left(r_{ij,L}^6\right)={\displaystyle \frac{s_6}{1+e^{-d(r_{ij}/\left(s_R{R}_{0ij}\right)-1)}}}.\end{array}$$

(5)

To keep things simple, only contributions (sum over $\mathbf{L}$) from neighboring cells were taken into account. The parameters for the damping function $s_6$, d, $s_R$ and $R_{0ij}$ were chosen according to the recommendations in the VASP manual [47].

Moreover, since the present work employed the TS vdW method, which essentially differs by the way how the dispersion coefficients $C_{6ij}$ are calculated (from the charge density instead of constants), the values were taken as calculated by VASP. Nonetheless, following the derivation in Ref. [24], it is possible to show that the variation of the dispersion coefficients is expected to be negligible. The dispersion coefficients within the TS method are calculated according to the following equations:

$$\begin{array}{c}\hfill {\alpha }_i={\nu }_i{a}_i^{\mathrm{free}},\end{array}$$

(6)

$$\begin{array}{c}\hfill C_{6ii}={\nu }_i^2{C}_{6ii}^{\mathrm{free}},\end{array}$$

(7)

$$\begin{array}{c}\hfill C_{6ij}={\displaystyle \frac{2C_{6ii}{C}_{6jj}}{\frac{{\alpha }_j}{{\alpha }_i}{C}_{6ii}+\frac{{\alpha }_i}{{\alpha }_j}{C}_{6jj}}},\end{array}$$

(8)

where ${\alpha }_i^{\mathrm{free}}$ is the polarizability of an isolated atom i and $C_{6ii}$ the dispersion coefficient for a pair of those atoms. Due to a assumption of a linear dependence between polarizablilty and effective volume [48] of the atom ${\nu }_i$ (obtained by employing Hirshfeld partitioning [49]), the effective volume is the only quantity which is calculated from the charge-density in order to obtain the dispersion coefficients. Therefore, plugging Equations (6) and (7) into Equation (8) reveals the linear dependence between ${\nu }_i$ and $C_{6ij}$

$$\begin{array}{c}\hfill C_{6ij}=\underset{\mathrm{constant}}{\underbrace{{\displaystyle \frac{2C_{6ii}^{\mathrm{free}}{C}_{6jj}^{\mathrm{free}}{\alpha }_i^{\mathrm{free}}{\alpha }_j^{\mathrm{free}}}{C_{6ii}^{\mathrm{free}}{\left({\alpha }_j^{\mathrm{free}}\right)}^2+C_{6jj}^{\mathrm{free}}{\left({\alpha }_i^{\mathrm{free}}\right)}^2}}}}{\nu }_i{\nu }_j.\end{array}$$

(9)

Here, ${\alpha }_i^{\mathrm{free}}$ and $C_{6ii}^{\mathrm{free}}$ are the polarizability and dispersion coefficient of an isolated atoms i, respectively, and can be found tabulated [50]. The dispersion coefficients for the POP-graphene differ only in a range of $\approx 1.7\%$$\frac{|max\left(C_{6ij}\right)-min\left(C_{6ij}\right)|}{〈C_{6ij}〉}$, caused by only small variations $\frac{|max(\nu )-min(\nu \left)\right|}{〈\nu 〉}\approx 0.85\%$ of the effective volume. Consequently the contribution arising from the changed electronic structure is rather small.

Comparing Figure 4b (calculated according to Equations (4) and (5) Figure 4c (full TS-VASP dispersion) suggests that the essential contributions to the adsorption surface topology are captured by the analytical form of Equation (4). Here, the geometry enters two terms: on the one hand directly through the $1/r^6$ term, and secondly through the damping coefficient. Setting the dispersion as well as damping coefficients to a constant value, Equation (4) simplifies to

$$E_{disp}\propto -\sum _{i=1}^N\sum _{j=1}^N\sum _{\mathbf{L}}{\displaystyle \frac{1}{r_{ij,\mathbf{L}}^6}}.$$

(10)

Thus, in Figure 4d the topology corresponding to only the geometry of the POP-graphene is shown, as given by Equation (10) (minimum-to-maximum scaled to fit the colour bar of the other dispersion maps). Clearly, the adsorption features such as preferred adsorption sites are related simply to the geometrical arrangement of the atoms involved in evaluation of the van-der-Waals dispersion, the major contribution to the adsorption energy.

5. Conclusions

After showing that our DFT-TS study is in good agreement with previous work [4], we could show that point defects such as single vacancies and Stone–Wales defects weaken the H${}_2$ adsorption in their closest neighbourhood, however, the effects remain spatially confined. On the contrary, the impact of a $(1,0)$ grain boundary on the adsorption properties is widespread, thus showing superior absolute values and large regions with enhanced adsorption properties as compared with the pristine graphene. POP-graphene, an example of a graphene derivative, outperforms even the grain boundary in terms of the enhanced area and the adsorption energies (adsorption is stronger). Furthermore, we showed explicitly by analyzing the effective atom volumes of POP-graphene, that the influence of the electronic structure to the dispersion correction has only a negligible impact. Hence no qualitative differences can be expected when employing the DFT-D3 method instead of the DFT-TS one, being also the reason for the aforementioned good agreement with previous literature on graphene. As a consequence, the changes in the adsorption capabilities mainly arise from the different geometrical arrangement of the carbon atoms in the investigated structures.