Abstract
Graphene-based nano-porous materials (GNM) are potentially useful for all those applications needing a large specific surface area (SSA), typical of the bidimensional graphene, yet realized in the bulk dimensionality. Such applications include for instance gas storage and sorting, catalysis and electrochemical energy storage. While a reasonable control of the structure is achieved in micro-porous materials by using nano-micro particles as templates, the controlled production or even characterization of GNMs with porosity strictly at the nano-scale still raises issues. These are usually produced using dispersion of nano-flakes as precursors resulting in little control on the final structure, which in turn reflects in problems in the structural model building for computer simulations. In this work, we describe a strategy to build models for these materials with predetermined structural properties (SSA, density, porosity), which exploits molecular dynamics simulations, Monte Carlo methods and machine learning algorithms. Our strategy is inspired by the real synthesis process: starting from randomly distributed flakes, we include defects, perforation, structure deformation and edge saturation on the fly, and, after structural refinement, we obtain realistic models, with given structural features. We find relationships between the structural characteristics and size distributions of the starting flake suspension and the final structure, which can give indications for more efficient synthesis routes. We subsequently give a full characterization of the models versus H2 adsorption, from which we extract quantitative relationship between the structural parameters and the gravimetric density. Our results quantitatively clarify the role of surfaces and edges relative amount in determining the H2 adsorption, and suggest strategies to overcome the inherent physical limitations of these materials as adsorbers. We implemented the model building and analysis procedures in software tools, freely available upon request.
Export citation and abstract BibTeX RIS
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
Since the rise of graphene [1], the problem of its production in forms useful for its employment in high-tech fields has received a great attention [2, 3]. In fact, its two-dimensional nature combined with the hexagonal symmetry of its lattice makes it a conductor with exceptional properties, useful for nano- and opto-electronics, but for other applications, its being flat and defect-less can be a drawback. For instance, for catalysis [4], gas sorting, water purification [5], the most relevant property of graphene is its large surface per unit mass (or specific surface area (SSA)), but porosity is needed as well; moreover, for applications in gas and energy storage [6] bulk properties such as the accessible volume per unit mass (or specific pore volume, SPV) are determinant parameters. Therefore for these applications, the large SSA must be conjugated to a three dimensional (3D) structure [7–9]. A possibility to realize these conditions are the graphene-based nano-porous materials (GNM) obtained using nano-flakes as precursors. These are quite a heterogeneous class including structures with a relatively high amount of sp2 carbon hybridization, characterized by different degrees of disordered porosity and in some cases by the presence of defects and adatoms, depending on the production method. At variance with their counterparts with porosity in the μm scale, obtained generally using nano-to-micro particles as templates [10], scaffolds with porosity strictly within the nano-scale are obtained directly from suspensions of nano-flakes exfoliated from graphene oxide [11, 12] possibly reduced [13], leading upon precipitation to extended networks of interconnected sp2 areas with SSA in the range of 650–1500 m2 g−1 [14–16]. SSA is then evaluated by the N2 adsorption on the basis of Brunauer–Emmett–Teller (BET) theory [17] (ideal graphene value is 2630 m2 g−1), and the chemical activation with KOH allows increasing it up to 3300 m2 g−1 [18, 19] and obtaining SPV up to 2.2 cm3 g−1 [20]. Based on a theoretical model consisting of perforated multilayers [19], this counter intuitive result was explained by the fact that the edges and perforations can adhere the gas, effectively increasing the SSA. In these materials, the best measured value of H2 adsorption is about 6% of gravimetric density (GD, the mass percentage of adsorbed H2) at 77 K and 20 bar [21].
In these disordered scaffolds, the structural measurement are limited to the pair distribution functions (PDF) [22] and average pore size (APS) [23] leaving quite a large amount of under-determinacy, and, consequently, high levels of arbitrariness in the model. The strategies chosen so far to produce models have been generally top-down or bottom-up. One example of the first class is the already mentioned 'perforated graphene' [19], using regular stacking of flat flakes. Other studies use either ideal [24–28] or defected multilayers [29], ordered 3D models as the open-carbon-frameworks [30], the carbon honeycomb [31] or the pillared graphene [32–34]. These models are useful to explain in first approximation the observed H2 adsorption capacity [35–37], but retain regularity and ideality to some extent. The lack of randomness in the structure prevents its realistic representation and the possibility of an accurate study of the adsorption, diffusion and interaction mechanisms between medium and adsorbate, which is of paramount importance to guide the design of devices for storage, sorting or catalysis [9, 18]. On the other hand, the bottom-up strategies based on molecular dynamics (MD) using carbon atoms as precursors can more easily incorporate disorder: starting from a random distribution of carbon atoms in gas phase, simulated annealing cycles [38] lead to porous scaffolds. Different morphologies can be obtained by changing the annealing conditions [39] (temperature, pressure [40] or density [41]), equivalent to changing the experimental conditions of production. This procedure relies on the use of force fields for the carbon interactions capable of describing the hybridization and bond order [42–45] and the final results depends on their accuracy. An alternative point of view is the reverse Monte Carlo method [46], where the atoms configurations are generated randomly and optimized until the simulated PDF matches the experimental one, in principle, independently from carbon interactions, although energetic or geometric restrains are often introduced [47] to avoid the formation of unrealistic structures and an excessive bias towards given experimental determinations. Bottom-up methods are more computationally expensive, preventing the generation of large models needed to mild the effect of boundary conditions, particularly critical in non-periodic systems, and a comprehensive throughput screening for the storage capacity, which needs a significant statistics of diverse models to cover a major number of cases of different SSA and SPV. Additionally, we remark that while pre-determining SSA and SPV is quite straightforward in regular models, this is not the case in disordered ones.
In this work we illustrate and implement a strategy that combines the advantages of top-down and bottom-up methods to the aim of efficiently building realistic models for GNMs with designed structural properties. The key ingredient is to follow an algorithm that resembles the experimental assembly process [48, 49]. Starting from precursors mimicking the exfoliated flakes allows sparing computational resources, while disorder is introduced both in the flakes distribution and in the subsequent optimization steps mimicking the treatments on the samples. We generate and analyze a statistical set of models spanning a large portion of the structural space, exploring the whole experimental range and beyond, and study the relationships between structural parameters. We then test the structures against their H2 adsorption capacity and derive quantitative rules to predict the hydrogen uptake with the aid of machine (ML) algorithms [50].
The next section reports the model building procedure, which is already a new result of this work, yet exploiting some standard simulation techniques. Section 3 reports the results on the structure characterization and the simulations of H2 adsorption. A summary and conclusions follow.
2. The model building procedure and other methods
The model building procedure consists of three steps, roughly corresponding the synthesis phases. The first is the preliminary model building by randomly mixing a distribution of flakes, corresponding to the exfoliation-dehydration-phase. In the second step the structure is optimized by means of MD, corresponding to environmental annealing. In the third step edges and other reactive sites are saturated with hetero-atoms (e.g. hydrogen), corresponding to chemical decoration or exposure of the structure to gases or other substances.
2.1. Step 1: preliminary model building
A scheme of the algorithm to build the preliminary model is reported in figure 1 (left) and consists of two modules: initialization (stages 1–5) and refinement (stages 6–7). In the initialization stage, flakes are generated with random size and placed with random location and orientation, first a set of seed of larger flakes (SF) and then a set of smaller filler flakes (FF). The shape and size are generated according to a user-defined distribution (defined by maximum and minimum sizes) that can be chosen e.g. based on those measured in the distribution in the flake suspension (see the SI is available online at stacks.iop.org/NANO/32/045704/mmedia for details). The SF placement defines the main morphology of the system (see figure 1(a), red), the FF (figure 1(b), in green) serve to give homogeneity and stability. Two parameters are additionally used to control the structure: d1 defines the minimum between the SFs and the FFs, and d2 defines the minimum distance between the centers of FFs.
Figure 1. Scheme of the algorithm to build the models. On the left, the step 1 (red), separated in two main modules: the initialization (1–5) including the mixing of flakes (inset (a) and (b)), and their perforation (c), and the refinement (6–7, (d), (e)). Step 2 on the upper right (blue) includes the structural optimization with two different force fields: AIREBO (f, f') and Tersoff (g, g'). Step 3 (right bottom, green) includes the H-decoration of edges, and optimization with MD (h, h').
Download figure:
Standard image High-resolution imageThe placement of flakes is iterated to adjust the density ρ to values slightly exceeding the target one, because in the following stages all the structural refinements are performed by deleting atoms. If ρ is lower, the procedure is repeated increasing the number of the initial SFs ('Check-Density' in figure 1). The steric hindrance at intersections among the flakes is eliminated (stage 4, figure 1 and the SI). The creation of perforations follows (figure 1(c), step 5) by deleting circular sections centered on casually chosen atoms. The size of perforation dp is randomly spanned between zero and a maximum value (currently set at 70% of the flake size) to avoid the scaffold collapse; to the same aim, during the process, restrains are used on the distance of the perforation from the edges (see SI.1.1). As a consequence, the average value of dp depends on the limiting values of the flake size, and either one or the other can be considered as inputs. This procedure generates structures whose value of perforation p (i.e. the ratio between the perforated and total area) normally ranges between 10% and 50%. Structures with p = 0 and p as large as 80% are also generated by forcing the perforation distribution to the upper or lower limits.
The refinement module basically adjusts ρ to the target value. Steps 6, and 7 ((d) and (e) in figure 1), consist in cleaning up the isolated atoms and too small flakes and in creating new perforations in sites at high local density. This cycle is repeated until the target ρ is reached. If ρ before the refinement is smaller than the target one, then the procedure is repeated from the SF placement. The outcome of this procedure is a model made of intersecting and interlacing perforated flakes randomly distributed in shape and orientation, with a given density and with controlled distances (figure 1(e)).
2.2. Step 2: structural optimization
In the second step the structure is endowed with physical character: the edges and intersections are relaxed, and the curvature of the flakes naturally forms. This is achieved by means of a multi-step protocol based on simulated tempering and annealing in which the temperature is alternatively raised to 1200 K and slowly decreased to 300 K during a total of several hundreds of ps (Simulation protocols reported in SI.1.2). MD simulations are performed treating the interactions with the AIREBO force field [43] (figure 1(f)) as implemented in LAMMPS [51], capable of dynamically generating a consistent network of bond orders with realistically hybridized sites and geometries [39, 52–55]. Flakes spontaneously bind, forming the scaffold geometry, free edges reconstruct and, where the local density is very low, structures as linear carbine chain can form (figure 1(f')). While these features realistically represent a carbon-only structure, to the aim of facilitating the subsequent step of decoration, non-reconstructed edges are more practical. Therefore, the structure is additionally subject to local optimization step with the Tersoff potential [42] (figures 1(g), (g')), which destabilizes the sp1 or distorted conformation in favor of the sp2, preparing the model for the decoration, without changing the global structure. The outcomes of this step are therefore two different models: one representing the carbon-only structures, and the second with more reactive edges to be used for the next decoration phase.
2.3. Step 3: decoration and functionalization
After a further clean-up of the structure to eliminate isolated carbon atoms or small clusters the structures are saturated with decorating atoms (figures and details in the SI.1.2). Currently, only H edge saturation is implemented (figures 1(h), (h')), being the most likely process occurring during H-storage processes, given the edges reactivity. However, partial saturations, decoration with other elements or even functionalization of other kind can be considered in next implementations, provided reliable FFs are available. The decorated scaffold is further relaxed with MD [41, 56].
2.4. Simulations setup and analysis methods
The structures are characterized using the mass densities (ρ, total and partial for hydrogen and carbon), the SSA, evaluated as in the BET theory, i.e. as the accessible surface to a probe of radius 1.7 Å normally used for the N2 molecule, and with a standard set of vdW radii and the SPV, both evaluated with the software Poreblazer [57], also used to evaluate the pore size distribution (PSD), while the radial distribution function g is evaluated with VMD [58], also used for visualization. An additional structural parameter related to perforation is the edge to surface ratio (ESR), measured as the relative amount of C atoms on the edges (NCedge/NC).
The H2 uptake is evaluated in the generated structures using Grand Canonical Monte Carlo (GCMC) at 77 K and 20 bar, considering scaffolds and the H2 molecules as rigid, and interacting with a Lennard Jones potential. GD is evaluated as the percentage of up-taken H2 (i.e. the edge decorating H is not counted as up-taken) with respect to the total mass, considering both the absolute and the excess uptake, that is
being Vpore the total pore volume, 
the H2 mass and 
its density evaluated in a void box with the same simulation parameters. Parameters definitions and relationships are reported in table 1. Additional details on the simulation protocols are reported in the SI.
Table 1. Structural parameters definitions and their relationships.
| Quantity | Definition | Evaluation method | Notes |
|---|---|---|---|
| ρ (g cm−3) | Total density including edge saturating atoms |
| |
| ρC (g cm−3) | Carbon-only density |
| Graphite: ρG = 2.26 gcm−3 |
| ρH (g cm−3) | H density (only edge saturating atoms) |
| |
(g cm−3) | H2 density (only gaseous part) |
| |
 (g cm−3) | Excess H2 density |
|
 = density of H2 reservoir at the same thermodynamic conditions |
| SSA (m2 g−1) | Specific surface area | The accessible surface to a probe having the N2 radius | Graphene: SSAg =2630 m2g−1 |
| SPV (cm3 g−1) | Specific pore volume | The accessible volume, defined by vdW spheres of carbon |
|
| δ = 'thickness' | |||
| APS, APV (Å, Å3) | Average pore size/volume | The average pore size obtained averaging the PSD and the corresponding volume | |
| PSD | Pore size distribution | Distribution of pore diameters | |
| ESR | Edge to surface ratio | Ratio between the edge carbon atoms (CH) and the total carbons | ESR = NCedge/NC |
| G(r), g(r) | Pair distribution function |
| G is the distance enhanced version of g |
| GD, GDex (%) | Gravimetric density |
|
|
The limited-memory Broyden–Fletcher–Goldfarb–Shannon machine learning algorithm [59] is used to relate the structural parameters to the H2 uptake, using 90% of the simulations as training set and 10% for validation. The neural network has been implemented with the MLPRegressor function included in the Python library Scikit-learn [60]. The mean absolute errors of our predictor are of 0.0039 over the training set and of 0.0044 over the validation set for the SPV, while of 0.039 over the training set and of 0.043 over the validation set for the hydrogen excess uptake (see SI.3).
3. Results: simulations and analyses
We generated 930 models of structures with hydrogenated edges and 205 with non-hydrogenated edges, obtained removing H from a subset of the former and relaxing. We first studied the relationships between the structural parameters. We then evaluated the H2 adsorption and its relationship to the structure.
3.1. Structures generation and analysis
Models were generated within a cubic supercell of 125 Å side, with periodic boundary conditions, using SFs up to 100 Å of size, with 11 different values of densities equally spaced in the range 0.2–0.7 g cm−3. We considered as upper limit the density of graphite, while the lower limit spontaneously emerges during the algorithm execution, as that below which instabilities and inconsistencies appear. Interestingly, the lower density lies in the typical range of metal oxide frameworks (MOFs, see figure S.13 in the SI and [16]). For each set of parameters two models were built to test the reproducibility of the algorithm. The main inputs are ρin, (d1 , d2) (hereafter globally named d), and p, which are, however, not all independent, and obtained with different combinations of accessory parameters such as average number and sizes of flakes and perforations. Table S.1 in the SI reports the inputs of the whole structure set, classified by ranges of the parameters.
Samples structures of the final models with hydrogenated edges are reported in figure 2 (a more complete set is the SI, figure S.4) for four representative values of the density, spanning the whole accessible range. The structures are characterized by pores and channels separated by sheets, or thin ribbons for lowest densities. Both the average value and the width of the PSD increase as ρ decreases (figure 2(e)) as a consequence of the relationships between ρ and the other structural parameters (see below), while at short distances (<1 nm) the PDF do not show large differences, being dominated by graphene-like features i.e. the peaks corresponding to the honeycomb symmetry, (black in panel (f)), broadened by the presence of disorder and of seven and five-member rings and structural defects. In fact, in good agreement with previous works [61], we find structures dominated by sp2 hybridization, with about 1% of sp3 at low ρ, increasing to 4% with ρ due to the larger amount of flake intersections (figure 3(c)).
Figure 2. Sample models of hydrogenated GNMs at different densities ((a)–(d), orange, red, green and blue label, respectively; densities reported in panel (e)/(f)). The SSA and SPV values are reported within the panels. Panels (e) and (e) report the pore size distributions and the radial distribution functions in corresponding colors. The black line in panel (f) corresponds to the PDF of graphene.
Download figure:
Standard image High-resolution image
Figure 3. (a) SSA versus SPV in the whole set of generated structures. Colors from red to blue are for increasing value of the flake distance and different values of the perforation (as reported). (b) Red, unit on the left: SPV-density plot. The solid line is the analytic excluded volume formula SPV = 1/ρ − (4π/3)(rC 3/mC), assuming all C atoms as isolated with rC = 1.7 Å the vdW carbon radius. The dotted line is a fit, leading to rC = 1.35 Å. Green: average pore size versus density. (c) The sp3 amount as a function of density and of edge to surface ratio. The dots are colored according to the density as indicated in the plot. Lines are fits with the formula sp3 = 0.017ρ2 – 0.03 ESR + 0.022. (d) SSA versus ESR, colored according to density (same color coding as c). The thin fitting lines are obtained with the formula in the text, with α = 2.5, A = 3.7 B = 0.45. The thick line is with the limit with ρ = 0.
Download figure:
Standard image High-resolution imageFigure 3 reports a quantitative analysis of the relationships among the structural parameters. Panel (a) reports the SSA versus SPV in the final structures. Accumulation lines are clearly visible, corresponding to the 11 different values of ρ. As the final density ρout is the result of subsequent adjustments, it can differ from the input ones ρin; the difference between the two is reported in the inset, and turns on average 1.5%, confirming the capability of our algorithm to control ρ. The lines can be understood considering that the accessible volume is roughly the total one diminished by the one occupied by the scaffold. For a system of isolated flakes this leads to the relationship
being δ the thickness of the sheets. In fact, the slopes of the lines are nicely fitted by the value δ = 3.4 Å that is, the vdW diameter of carbon, a reasonable approximation of the thickness (fitting lines reported in black in figure 3(a) for a few values of the density, as indicated). In principle a density ρfit can be fitted from the intercepts. This coincides with ρin and ρout at low densities, while it is larger at high densities where the sheets intersections are no more negligible.
At each fixed ρ, the range of possible morphologies is spanned by varying the average flake distance d (d1, d2) and the average perforation p in order to span as much as possible the SSA and SPV values. The lower limits (red dots) are obtained with low value of p and large value of d. By decreasing d, one can climb up to the upper limit (blue dots). Accordingly, the value of p changes to maintain the input ρ. It turns out that the lower range limit corresponds to scaffolds with large, rather distant and almost defect-less sp2 areas and consequently lower SSA, while the upper limiting structures are frameworks made of thin ribbons produced by a large number of defects, similar to the MOFs, with smaller d and larger p values to maintain the scaffold stability.
Figure 3(b) reports the dependence of the porosity on the density of the final structures. The SPV (in red) decreases as the density increases, as effect of the increase of the volume fraction occupied by the atoms. This is always slightly larger than the free volume evaluated considering isolated spheres with the carbon vdW radius (the red solid line is with rC = 1.7 Å) due to the atoms superposition and the flakes intersections. A better fit is obtained using a smaller effective radius (dotted line, corresponding to rC = 1.35 Å). At any given density, therefore, the SPV varies little, while the APS can span a relatively large interval from 5 Å at large density, to 20 Å for small densities (figure 3(b), green dots). In fact, as the density decrease, the sheet size, distance and perforation size can increase.
The dependence of sp3 amount on the ESR and ρ (color coded) is reported in figure 3(c), and is seen to increase with ρ due to the relative increase of intersections and decrease with ESR since the presence of perforations reduce the intersection possibility. As shown by the fitting lines, the dependence on ESR is linear in first approximation, while the increase with ρ is quadratic (fitting parameters reported in the caption). The fit is better at lower densities.
Finally, the dependence of SSA on ESR is reported in figure 3(d). The dependence is fitted by a bundle of lines whose slope decreases as the density increase, reported as thin black lines in the plot. This behavior can be explained as follows: if the amount of intersections is small (i.e. for low ρ) the SSA can be approximated by that of graphene augmented by the presence of edges, i.e. linearly increasing with ESR. This explain the low density behavior of SSA (colors magenta to yellow) and is consistent with previous works on perforated graphene [19]. The intersections however, acts the way round, decreasing the effective SSA. Overall, we found that the whole set of data is nicely fitted by a three parameters function
where the parameter α = 2.5 is the enhancement contributed by the edges to the SSA, while A and B describe the decrease due to the intersections (numerical values of A and B in the caption). The low density behavior of SSA enhancement is obtained with ρ → 0. At high density the intersection effect prevails and supersedes the edge effect, and overall the SSA decreases. The whole analysis is valid both for the systems with hydrogenated edges (reported in figure 3) and for the non-hydrogenated-edges systems, reported and compared in the SI (figure S.2). The fitting functions are the same, with very similar parameters, one remarkable exception being the parameter α which assumes the larger value 2.7, consistent with the larger adsorption capability of carbon reconstructed edges with respect to the hydrogenated ones. The full set of input/output data are available upon request.
3.2. Hydrogen adsorption
The hydrogen adsorption was evaluated in all the generated structures. Samples of loaded structures taken along the upper and lower boundaries of explored regions in SSA-SPV plane are reported in figure 4. Panel (a) reports one high loading case, taken in the 'MOF-like' region, in the upper right part of the SSA-SPV plane, with low density and high ESR (values given in the caption). As it can be seen from the snapshot, H2 (in orange) adheres to ribbons and edges, and fills nano-sized concavities. Panel (b) is a case of rather large ESR but low high density. Due to the occupation of the volume by the carbon, the SSA is much reduced and the GD is halved. In (c) we report an opposite case, of low density (as in (a)), but also low ESR, i.e. large sheets with low value of perforation, which makes these structures 'foam like' and lye in the lower edge of our explored region in the SSA-SPV plane. In this case the adsorption is still reduced with respect to (a) due to the absence of edges and assumes an intermediate value between (a) and (b). In (c) we adopted a different visualization to highlight the presence of layers: from the whole GCMC trajectory, we generated two volume density maps for H2: one for the 'first layer' defined by all hydrogen molecules nearer to C atoms less than 4 Å, and the other for those in the interval 4–8 Å from carbon, defining the 'second layer'. The first and second layer maps are visualized as iso-density surfaces, in yellow and red, respectively. As it is apparent, while the first layer is clearly visible and well formed, the second layer is sporadic and fragmented, and present mainly in the concavities or elbow like regions. This tendency of H2 to cumulate in concavities was previously observed [62], and attributed to an enhancement of vdW interactions. A more exhaustive set of H2 loaded structures is reported in the SI, figures S.9–S.11.
Figure 4. Scaffolds (hydrogenated-edges series) loaded with H2. (a) ρ = 0.2 g cm−3, ESR ∼ 0.5, (GD ∼ 4.7%). Scaffold atoms in grey, H2 molecules in orange. (b) ESR ∼ 0.4, ρ = 0.7 g cm−3, (GD ∼ 2.5%). Same visualization as in (a). (c): ρ = 0.2 g cm−3, ESR ∼ 0.25, (GD ∼ 3.5%). Hydrogen molecules are represented by means of iso-density surfaces of H2 molecules in the first layer (yellow) and in the second layer (red). Un-layered hydrogen is colored in mauve.
Download figure:
Standard image High-resolution imageWe quantitatively explored the relationships between adsorption and the structural parameters (figures 5 and 6). We make first two general observations: first, GDex is always less than GD and the difference between the two increases as ρ decreases, as an effect of the increase of void space. Second, GD (GDex) for the hydrogenated-edges systems is less than those for the non-hydrogenated-edges ones, the former spanning the interval 1.7%–7% (1.4–4.5) the latter the interval 3–7.6 (2.3–5.2). Again, this is consistent with the larger adsorption capability of carbon reconstructed edges with respect to the hydrogenated ones. These results are compatible with previous theoretical and experimental data [19] (see also figure S.13 in the SI).
Figure 5. Dependence of H2 adsorption on the structural parameters for the edge hydrogenated scaffolds. (a)–(c): GDex versus density, APS and ESR, respectively, colored according to SSA (color coding reported in (a)). (d), (e) GDex versus SSA and SPV, colored according to ESR (color coding reported in (d)). (f) and (g), SSA versus SPV and ESR, respectively, colored according to GDex (color coding reported in (f)). The lines in (g) are eye-guides to separate the areas with increasing GDex. Their slope is 4000 m2 g−1 and the intercepts are evenly separated by 380 m2 g−1 between 900 and 2420 m2 g−1. (h) Scatter plot between the simulated GDex (abscissa) and the estimated Gex (ordinate) obtained with the formula GDex (%) = 0.0013 (g m−2) SSA – 5.26 ESR + 0.85.
Download figure:
Standard image High-resolution image
Figure 6. Same as figure 5, for non-hydrogenated-edges scaffolds. Same color coding and symbols. In (g) the lines slopes is 3000 m2 g−1 and the intercepts are separated by 450 starting from 1200. In (h) the formula used for GDex estimate is GDex (%) = 0.0010 (g m−2) SSA – 3.0 ESR + 1.5.
Download figure:
Standard image High-resolution imageIn fact, figures 5 and 6 clearly show an inverse dependence from density (panels (a)) and a direct dependence on APS, SPV and especially SSA, showing a linear correlation and therefore appearing the main determinant of GDex (panels (b), (e) and (d) respectively). However, especially in the case of the hydrogenated edges (figure 5(d)), the linear dependence appears to be strongly modulated by ESR, as shown by data coloring. This effect is better seen in panels (g) of figures 5 and 6, where the scatter plot SSA versus ESR is colored according to GDex: different uptakes areas are clearly separated in parallel strips, indicating that the uptake increases linearly along a direction roughly orthogonal to those strips. Therefore it is possible in first approximation to write an empirical relationship for the uptake depending on SSA and ESR as 
We fitted the parameters A, B and C on the data in panels (g), and then used the formula to predict GDex as a function of SSA and ESR. The comparison predicted versus simulated is reported in panels (h) of figures 5 and 6 and displays a relative root mean squared deviation of ∼1% and ∼2%, respectively, indicating that these linear combinations of the parameters SSA and ESR are very good determinants of the GDex. In both cases the main determinant is SSA, but, while in the non-hydrogenated edges cases the contribution of ESR can brings a correction always less than 1% in GDex, in the case of the hydrogenated edges this contribution is larger, sometimes exceeding 2%.
While the direct linear dependence of GD on SSA is not surprising, the negative correction by ESR, measuring the amount of edges, is counter intuitive, as it was previously shown in the literature that the presence of edges increases hydrogen adsorption. We additionally show in figures 5, 6 panels (c) the GDex generally increases with ESR, except at high densities (blue color) when intersections become predominant. These apparently contradictory observations can be reconciled considering the definition of SSA adopted in this work (and generally in the theoretical works): SSA is evaluated as the accessible surface to a probe with the vdW parameters of N2. This method does not discern where molecules are adsorbed, or the differential adsorption capability of edges, especially if decorated. Edges are assumed to adsorb as the surfaces, and bring a relative enhancement to SSA only due to the larger exposure of edge atoms with respect to surface ones. This turns out an overestimate, especially in the case of H-edge saturation, since the adsorption capability of H is less than that of C. Therefore, a negative correction proportional to the edge amounts emerges. In fact, the correction is smaller in the case of reconstructed C edges, though still present, indicating that, overall, even the adsorption capability of C on edges is not as larger as it should be as an effect of the larger exposure.
Two important remarks naturally stem from these results: the definition of SSA generally used in theoretical works differs by that used in experimental determinations, since the latter is strictly proportional to GD and therefore already includes the differential adsorption of edges, which, as seen, is present also in carbon-only structures, though smaller than in the hydrogenated-edge case. This also imply that simplified theoretical approaches based only on the geometrical determination of SSA can bring systematic errors in the GD evaluation. On the other way round, it also clearly emerges that the decoration of edges with highly adsorbing elements could bring a correction of opposite sign, i.e. increase the adsorption, and this could be used to design scaffold with tailored GD.
3.3. The effect of edge hydrogenation
While in the previous section we considered two extremal cases, i.e. completely hydrogenated edges and completely dehydrogenated ones, the real systems are more likely to undergo partial edge hydrogenation, during the first phases of exposure to H2, even starting from carbon-only structures. To evaluate the effect of a partial hydrogenation on H2 physisorption we used a heuristic approach completely based on machine learning. We used as training set the values of GDex in all structures. The inspection of prediction errors (see SI.3) let us single out the most significant input variables among the different correlated descriptors of our system, leading to an engineered function of SSA, total density and the ratio of number of edge chemisorbed H over the whole number of carbon atoms, hereafter referred to as H/C.
Based on that, we are able to predict the value of observables of interest as continuous functions of selected input variables within an interval larger than the simulations data range (see SI section S.3). In figure 7 we show the values for GDex predicted by the neural network in the continuous SPV-SSA range and for a set of given intervals of H/C. The function correctly reproduces the training and validation data (superimposed as dots), and produces interpolated and extrapolated data (even extending in the region of structural instability of the scaffold). The general trends confirm SSA as main determinant of H2 physisorption, but it is also confirmed the strong effect of edge hydrogenation, which should be minimized to improve GDex.
Figure 7. Prediction of the GDex as a function of SPV and SSA for different values of H/C (indicated over the plots), according to the ML algorithm. The shades of color describe the H2 uptake as specified in the legend. Superimposed dots correspond to training and validation data within the given H/C range (around 90% of them were used as training data, with the remaining part contributing to validation).
Download figure:
Standard image High-resolution image
4. Summary and conclusions
In summary, in this work we first design and implement an algorithm to build graphene-based scaffold with porosity in the nano-scale. We demonstrate that our method is capable of building realistic models with pre-determined density and porosity. This goal is achieved using a multi-step strategy mimicking the real synthesis. After creating a statistical set of edge hydrogenated and carbon-only structures with reconstructed edges, we analyze them and highlight the correlations between structural parameters (SSA, SPV, ρ and ESR) and which could be used to interpret the structural data on these somehow elusive systems.
We subsequently simulate the H2 physisorption in the whole structures set and study the structure-adsorption relationship, deriving quantitative empirical relationships between the structural parameters and GDex, which are able to predict the simulated adsorption with the precision of a few %. These relationships additionally identify as main structural determinant, beside SSA, also ESR, measuring the relative edges amount. Remarkably, we found that the ESR can bring substantial correction to the GD, as an effect of the differential adsorption of edges. We further investigate the effect of edges decoration with hydrogen by mean of a machine learning algorithm, which confirmed that the GD is indeed modulated by the amount of chemisorbed H on edges, acting in a negative way, i.e. decreasing the GD due to the reduced physisorption capability of hydrogen.
Our results also indicate that using simplified models with non-reconstructed or non-functionalized edges and approximate determinations of structural parameters can bring non-negligible systematic errors in the adsorption evaluation. These, in turn, can be prevented by accounting for the differential adsorption of edges, depending on their structure and/or decoration. Additionally, our quantitative analysis shows that the differential adsorption of edges could be exploited to overcome the graphene the scaffolds limitations and design materials with tailored or selective gas physisorption, which could be obtained by decorating the edges with elements or functional groups specifically designed to be affine to H2 or other given gases.
Acknowledgments
This research was funded by EU-H2020, Graphene-Core1 (agreement No. 696656), Core2 (agreement No. 785219) and FETPROACT LESGO (Agreement No. 952068), by the Italian Ministry of University and Research (project MONSTRE-2D PRIN2017 KFMJ8E), by CINECA awards IsC61_MGchpDFA (2018) IsC69_EFaRe (2019), IsB19_DiNaGra (2019). We gratefully thank CINECA staff for technical support.