Adsorption-Induced Expansion of Graphene Oxide Frameworks with Covalently Bonded Benzene-1,4-diboronic Acid: Numerical Studies.

Graphene oxide frameworks (GOFs) are interesting adsorbent materials with well-defined slit-shaped pores of almost monodisperse separation of ∼1 nm between graphene-like layers; however, the exact nature of the structure has remained undetermined. Recently, GOFs were observed to swell monotonically upon the adsorption of methane and xenon under supercritical conditions. Here, we present the results of molecular dynamics simulations of the adsorption of methane and xenon for various proposed GOF structures based upon force fields based on ab initio B3LYP density functional theory calculations. The simulations reproduce well both the adsorption isotherms and the expansion of the interlayer spacing for methane and xenon for a model of GOFs formed by covalently bonded benzene-1,4-diboronic acid oriented at quasirandom angles with respect to the graphene layers.

Introduction

Carbon-based porous materials have been studied extensively for their potential for energy storage in the form of hydrogen or methane, other gas storage applications, and gas separation (e.g., sequestration of carbon dioxide). The reason for this interest in carbon nanostructures is their tremendous variety, which includes various amorphous activated carbons,[1−3] metal organic frameworks (MOFs),[4,5] covalent organic frameworks (COFs), and more recently graphene oxide frameworks (GOFs). In general, the storage capacity of these materials is roughly proportional to their specific surface area (SSA, which often is close to ∼2600 m2/g of graphene) and is maximal for pore sizes in the neighborhood of ∼1 nm. For maximal storage delivery, Bhatia and Myers[6] demonstrated that there is an optimal adsorption enthalpy of ∼18 kJ/mol for room-temperature operation, and Kuchta et al. demonstrated that an optimum adsorbent should be homogeneous.[7]

Graphene oxide frameworks (GOFs) are a particularly promising structure for various physicochemical processes. GOFs are formed when graphene oxides are linked by different functional groups,[8] thus producing a layered structure[9−11] with possibly tunable pore widths depending upon the type of binder used, which may be optimized, e.g., for storage of a particular gas.[9,12−15] Notably, similar layered materials have been produced by silylation reactions.[16,17] A detailed understanding of the structure of GOFs at the atomic level is thus desired to produce future high-performance nanoporous materials for applications such as gas storage, separation, catalysis, supercapacitors, batteries, etc.

In this paper, we focus on the structural properties of a type of GOF formed using an intercalation procedure originally proposed by Burress et al.[9] and reproduced by Mercier et al.[11] and Schaeperkoetter et al.:[18] the dissolution of benzene-1,4-diboronic acid (DBA) in methanol and then intercalating the DBA into graphite oxide (GO), which lead to GOFs with SSA of ≈400 to 1000 m2/g and interlayer distances of ∼9 Å.

Originally, there were three competing models of GOFs (possibly) linked by DBA molecules. Burress et al.[9] assumed that GOFs were composed of graphene oxide layers separated by rigid DBA linkers with oxygen atoms covalently bonded to both graphene planes and benzene diboronic linkers; we refer to this structure as (i) “covalent pillared GOF” (CP-GOF). Mercier et al.,[11] however, posited that the proposed structures would be unlikely given the characteristic size of the DBA linkers (i.e., the pores would be larger than observed experimentally), and that their observations of swelling of the pores by as much as 50% in methanol seemed unlikely given the strength of typical covalent bonds.[11] Thus, Mercier et al.[11] proposed two alternative models for the GOF-DBA structure: (ii) the “van der Waals GOF” (vdW-GOF), where covalent bonds are formed alternatively on alternating walls of the pore, and van der Waals interactions between the DBA tails securing the structure together and a (iii) “fluid-DBA” model (fDBA-GOF), where DBA molecules possibly interacting via enhanced electrostatic interactions with polar groups in the GO sustain the two GO planes apart. As we shall see, our analysis rules out models (i)–(iii) and proposes an alternative model: (iv) “covalent angled GOF” (CA-GOF), which is compatible with experimental data. The models are depicted in Figure .

Figure 1

Postulated models of GOFs: the covalent pillared GOF (top left), the fluid DBA (top right), the vdW-GOF (bottom left), and the covalent angled GOF (bottom right). A stabilizing layer is included at the bottom to simulate an underlying solid structure (see also refs[9−11, 18][9−11,18]

A significant breakthrough came in 2019 when Schaeperkoetter et al.[18] probed these structures by in situ neutron diffraction during supercritical adsorption of various gases at pressures in the 0–140 bar range. Thermal neutrons from the 10 MW University of Missouri Research Reactor were particularly well suited for structural investigations of GOFs[19] since they can readily penetrate the thick high-pressure aluminum sample cell (in contrast, X-ray scattering cross section increases with Z2, making it very difficult to probe carbon structures surrounded by thick metal walls).

Schaeperkoetter observed that the interlayer spacing d001 increases gradually and monotonically as a function of pressure for methane, hydrogen, and xenon during supercritical adsorption.[18] The gradual change did not fit the typical either “gate-opening” transitions or other contraction/swelling of nanoporous materials that had been previously observed under subcritical adsorption[20−26]—in fact, in gate-opening systems, a bulk phase change usually precedes the gate-opening such as capillary condensation, but our gas’s isotherms were measured at supercritical conditions making a bulk phase change impossible. Another notable observation was that this swelling could be mapped onto a common expansion curve for all three gases using a scaling directly related to the van der Waals gas–solid interaction energy[18] in a manner similar to the law of corresponding states for gases.[27,28]

As we shall see below, our simulations of the adsorption of methane and xenon in the four GOF-DBA models proposed are consistent only with the CA-GOF model shown in Figure (bottom right panel).

Methodology

In this work, we analyze by means of molecular dynamics (MD) simulations the behavior of various models (Figure ) of GOF-DBA under conditions similar to those of the experiments of Schaeperkoetter et al.:[18] methane and xenon pressures in the 0–120 bar range at room temperature (T = 300 K). All simulations were completed using the NAMD2[29] MD code and analyzed using visual molecular dynamics (VMD) and various TCL scripts.[30] Interactions were cutoff for distances larger than 12 Å, and a neighbors’ list with a radius of 14 Å was maintained for accelerated calculations of forces. For each run, 5000 steps of energy minimization were performed, followed by 3,000,000 × 1 fs steps (i.e., 3 ns) of MD simulations in the canonical (N, V, T) ensemble using a velocity rescaling thermostat applied to all atoms (except for a single fixed stabilizing layer) every 20 fs with configurations saved for analysis every 1000 time steps (i.e., every 1 ps). The last 2 ns of each simulation was used for calculations of time averages (all simulations equilibrated in less than 0.5 ns). All simulations were performed inside a 21.3 Å × 37.2 Å × 50 Å parallelepiped box with periodic boundary conditions (PBC) in all directions; see Figure .

Figure 2

Head on view of a typical system along the x-axis (top left) showing regions (small regions (small and large y)) where gas molecules can transfer between the pore and gas phases. The system viewed along the y-axis (top right) shows carbon atoms in GOF bonded across the PBC. The “Gas Regions” in yellow, situated farther than the potential cutoff distance (12 Å) from the GOF, depict where the gas density is measured[32] and the pressure calculated using NIST’s Thermophysical Properties of Fluid Systems database.[31] The bottom panel illustrates the “strip” shape of the simulated GOF, i.e., infinite in the x- but finite in the y-direction.[32]

In each simulation, the GOF structure (models i, ii, iv), GO, and DBA (model iii) were simulated fully atomistically, and the atoms were allowed to move with the exception of a fixed graphene bottom layer that simulates a solid graphite substrate, increases the stability of the simulation, and facilitates the analysis by keeping the structures roughly aligned parallel to the xy plane. This is a reasonable model since GOF samples have SSAs of ≈400 to 500 m2/g[18] compared to graphene’s theoretical surface area of 2600 m2/g.

The gases (methane, xenon) were also simulated fully atomistically. The number of gas molecules Nmet ∈ [0, 400], NXe ∈ [0, 400] was varied, and the pressure calculated using NIST’s Thermophysical Properties of Fluid Systems database[31] from gas densities averaged far from the GOF (and DBA), i.e., in the bottom 10 Å of the simulation box; see Figure and ref (32).[32]

For interactions, we used the CHARMM22 force field[33]

The bond term is modeled in the harmonic formwhere b is the bond length between a pair (1–2) of atoms, b0 is the equilibrium bond length, and kb is the bond force constant. To avoid unnecessary short time scales, C–H bonds were made rigid (b = b0) using the RATTLE algorithm,[34] which allowed a time step of 1 fs. The second term accounts for the angle deformationwhere θ is the angle between a sequence of three bonded atoms, θ0 is the bond-angle equilibrium value, kθ is the angular harmonic stiffness constant, S0 is the equilibrium distance between 1 and 3 pairs, and kUB is the Urey–Bradley constant. The third term accounts for four-body dihedral torsion contributionswhere χ is the dihedral angle formed by a sequence of four bonded atoms, kχ is the torsional stiffness, n is the multiplicity factor, and δ is the phase shift. The final term makes up the nonbonded interactions derived from two-body interactions for atom pairs either in the same or a different molecule and adatom–substrate pair interactionswhich are applied only to atom pairs separated by at least three bonds,[35] with 1–4 interactions modified by a scaling factor of 0.4.[35,36] Here, q are the Mulliken partial atomic charges, r are the interatom distances, ε0 is the dielectric constant, ε is the Lennard-Jones potential depth, and rmin is the distance of the Lennard-Jones minimum. For heterogeneous atom pairs, we use the Lorentz–Berthelot combination rules[37]All LJ potentials were taken out to a pair separation of r = 10 Å and then smoothly diminished to zero at a cutoff distance rc = 12 Å. Coulomb interactions were calculated using the particle mesh Ewald (PME) summation method.[38]

In this work, the values of potential parameters are the result of a comprehensive study combining ab initio density functional theory (DFT/B3LYP)[39−42] calculations using the Gaussian09 code[43] (including a base study, geometry optimization, variation of geometric parameters to determine elastic constants) and the CGenFF database[44,45] for the dispersion terms. The methodology of the ab initio study and values of the interaction parameters used in this work are listed in the Supporting Information (see Tables S1 and S2; Figures S1 and S2).

Results

Model Descriptions

MD simulations of the adsorption of methane and xenon were performed for the models of GOFs described above (see Section and Figure ):

For model (i) CP-GOF9, a number of DBA molecules (3 or 4 in our simulation box, as per the experimental B/C ratios of GOF samples[9−11,18]) are covalently linked to both sides of the pore through C–O bonds (eq ) to graphene at sites that are almost exactly on top of each other (Figure , top left). This results in a layer separation d001 = 10.6 Å in disagreement with the experimental neutron diffraction results for d001: between 9.30 Å (in vacuo) and 9.75 Å (PXe = 40 bar)[18] (note that in ref (18)d001 = 10.3 Å using X-ray diffraction, but this result is for samples in “wet” air, consistent with observations in ref (11)). The covalent pillared GOF structures were unsurprisingly very stable and rigid and, as we shall see, showed no expansion during gas adsorption due to the substantial rigidity of the C–O covalent bonds and within the DBA itself against stretching.

For model (ii) vdW-GOF11, we also placed 2–4 DBA molecules (again consistent with experimental B/C ratios[9−11,18]) but bonded them through C–O bonds to alternate sides, allowing nonbonded (eq ) interactions between the DBA molecules to stabilize the structure (Figure , bottom left). For 2–3 DBA molecules, the structures were unstable (i.e., the pores collapsed to d001 < 6 Å, thus making them inaccessible to gas adsorption, see Figure , left panel), requiring at least 4 DBA molecules to remain open, but in that case, the d001 ≈ 11–12.5 Å was also too large (and, as we shall see later, remarkably pressure independent, as expected for larger pores[24−26]). To test if the problem was a finite-size effect, we doubled the xy dimensions of the simulation cell (quadrupling the GOF surface area), but the results were unchanged: the van der Waals GOFs resulted in either a collapsed pore or a very large and stable pore.

Figure 3

Typical collapsed van der Waals GOF (left) and fluid-DBA (right) models. These configurations almost completely excluded any gas molecules from adsorbing inside the pores.

For model (iii) fDBA-GOF, we created GO by adding epoxy and hydroxyl groups randomly to graphene based on the Lerf–Klinowski model[46−48] and randomly placed 2–6 DBA molecules between the GO layers (see Figure , top right). We tested numerous randomly produced GO/fluid-DBA configurations, but the results were always the same: a complete collapse of the pores to interlayer distances d001 < 7.8 Å that almost completely excluded any gas molecules from adsorbing inside the pore; see Figure , right panel.

Here, we propose a new fourth model (iv) CA-GOF, which retains the main characteristics of Burress’[9] but that is compatible with experimental observations of base layer spacing d001 ≈ 9.3 Å and, as we shall see below, shows adsorption isotherms and swelling of the interlayer spacing consistent with the experimental results of Schaeperkoetter et al.[18] for both methane and xenon. The generalization is intuitive: instead of linking the DBAs perpendicularly to the graphene planes, we first covalently bond them (eq ) to a random pair of C atoms in one layer of graphene (i.e., forming C–O bonds, it should be noted that the algorithm does not allow C–O bonds for different DBA molecules within the same graphene unit cell). Second, a random search is performed for a C–C bond in the second layer of graphene within the geometrical constraints to produce the desired layer spacings (i.e., the binding sites are displaced between 4.5 and 7.5 Å in the graphene plane); see Figure , bottom right. In a typical simulation, we place three DBA molecules per simulation cell so that the boron content is comparable to that observed.[9,11,18]

Initial energy minimization in NAMD2 is then performed to eliminate geometrically unstable configurations, in addition to eliminating those for which the starting pore spacing d001 (N = 0) differed by more than ±0.5 Å from the experimental results (9.3 Å). The idea is that the configurations that match the pore spacing of Schaeperkoetter et al.’s experiments[18] are representative of the GOFs from the said experiments.

Determination of Interlayer Spacing d001 and Variation with Methane and Xenon Adsorption

We created VMD/TCL scripts to calculate the average pore spacing by finding the density of graphene carbon atoms at specific z-values of separation between pores P(z). In the script, the first 1,000,000 MD steps of “thermalization” are discarded, and statistics was collected for the last 2,000,000 steps (2000 frames separated by 1 ps), as described in Section . Figure shows variations of P(z) distribution for some values of methane pressure Pmet for one particular configuration of model (iv) “covalent angled GOF.” We observe that the density distribution P(z) varies with the number of gas molecules (pressure) and, in general, we observe that the average pore spacing d001 = ⟨z⟩ increases with pressure as observed via neutron diffraction by Schaeperkoetter et al.[18]

Figure 4

Height density plot of carbon atoms in the top layer of graphene in the simulated GOF relative to the bottom layer at various methane pressures for one typical “covalent angled GOF” configuration. We observe a gradual increase of the average pore spacing d001 = ⟨z⟩ with pressure similar to what is observed experimentally.[18]

We should note that since our structures were created at random, some resulted in pores that partially collapsed (e.g., because all DBA linkers were clustered in one small area)—these runs, as mentioned earlier, were discarded—it is thus important to remark that any conclusion must be drawn from analyzing averages of large ensembles of configurations. Figure shows averages of d001 variations with methane pressure for the four GOF-DBA models discussed in this paper. It is evident that the interlayer distance d001 for CP-GOF9 and vdW-GOF11 are both too large (∼10.3, 11–12.5 Å), and for fDBA-GOF11 too small (7–8 Å, and excluding any gas from its interior) when compared to Schaeperkoetter’s neutron diffraction results (9.3–9.75 Å).[18] Only the proposed CA-GOF is consistent with experiments even before the small variation (in this scale) of d001 with methane pressure is discussed below.

Figure 5

Pore spacing d001 calculated by averaging over numerous configurations for each of the discussed GOF-DBA models (Figure ) and comparison with experimental results of Schaeperkoetter et al.[18] Error bars of the experimental data are omitted because they are too small at this scale. The dispersion of the computational data is shown in Figure .

Figure 6

Pore spacing relative increase d/d0 as a function of methane pressure Pmet. Thin lines are individual configurations color coded by their low-pressure variation in d001 vs the experiment[18] (note: d0 = 9.35 Å). Thick lines represent averages within each class.

In the case of the CA-GOF model proposed in this paper, the random arrangement of DBA molecules causes a range of possible distances between graphene sheets. For the distance between graphene sheets to be equal to the experimental distance, diborane bonds with two graphene layers must be shifted in the xy direction by ∼6 Å. For our analysis, we picked structures with random shifts between 4.5 and 7.5 Å and discarded results in which geometries resulted in pores that collapsed. Figure shows how the layer spacing d001 = ⟨z⟩ (see Figure ) varies with the methane pressure for 32 random configurations color coded by how close they are to the low-pressure experimental spacing.[18] It is evident that the results are quite dispersed. However, averaging over these configurations (thick lines in Figure ) results in curves that are in reasonable agreement with neutron diffraction experiments,[18] and unsurprisingly, the average calculated for the seven configurations closest to the experimental value has a better agreement than the more dispersed data set. A similar analysis was performed for the swelling of the interlayer spacing upon adsorption of xenon. The data set is more limited but shares the methane’s characteristics.

Figure shows the variation of the layer spacing d001 as a function of pressure during adsorption of methane and xenon. The lines represent the ensemble averages of many GOF configurations as described above, and symbols are from Schaeperkoetter et al.’s in situ neutron diffraction experiments.[18] Overall, the ensemble average variation of layer spacing with gas pressure for this model is in very reasonable agreement with experiments. The larger variation of d001 with pressure for Xe was explained in ref (18) from the different gas interaction parameters.

Figure 7

Pore spacing absolute increase of d001 as a function of methane and xenon pressure. Experimental data from Schaeperkoetter et al.[18] The larger variation of d001 with the pressure for Xe was explained in ref (18) from the different gas interaction parameters.

Adsorption Isotherms

We also created VMD/TCL scripts to compute the average number of gas molecules residing inside the pores to calculate the absolute adsorption isotherms for CA-GOF model structures (last 2000 frames, as in Section ). The number of adsorbed molecules inside the GOF structure is averaged, while separately the gas pressure is calculated using NIST’s Thermophysical Properties of Fluid Systems database[31] from gas densities averaged far from the GOF (and DBA), i.e., in the bottom 10 Å of the simulation box; see Figure (32)Figure shows the individual configurations’ adsorption isotherms and their ensemble averages (see discussion in Section , Figure ). Figure shows a comparison between the experimental and average simulated adsorption isotherms for methane and xenon. Xenon adsorbed noticeably more rapidly in the simulations than in Schaeperkoetter et al.’s[18] experiments, but overall there is a reasonable agreement for both gases considering the large degree of variability in the simulated GOF structures.

Figure 8

Absolute adsorption of methane in “covalent angled GOF” systems and their ensemble averages and experimental results.[18]Figure shows the color scheme.

Figure 9

Absolute adsorption of methane and xenon in “covalent angled GOF” ensemble averages for the ±0.1 Å tolerance level compared to experimental results.[18]

Discussion

Our observations of monotonic and gradual swelling of the d001 spacing for the CA-GOF model are consistent with Schaeperkoetter et al.’s observations.[18] Perhaps the main criticism of any of the covalently bonded models (model i or iv, see Section ) is by Mercier et al.:[11] they posit that because the pore size swells up to 15.4 Å when in water, the DBA cannot be covalently linked as 12 Å would be the limit of the height.[11] However, an alternative explanation for the large swelling when flooded with a polar solvent may be boronate esterification, a well-studied reversible reaction[49,50] where a boron bonded to two oxygen and a diol, in this case carbon, will have its oxygens detach from a surface such as graphene in the presence of water, leaving behind OH on the surface as the boron’s two newly detached oxygens each gain a hydrogen, all within a short time.[50] When the boronate’s oxygens bond back with the surface, the reaction produces water. If no water is present, the boronate remains bonded to the surface. Therefore, flooding the pores with water could, in theory, detach all DBA molecules and then swell the pores to distances too great for the DBAs to bond to both sides. Conversely, continuous heat-induced evaporation of water in a GOF will result in DBAs being only covalently bonded to the surface. In conclusion, it is possible for the covalently bonded models to be consistent with pore sizes swelling up to 15.4 Å.

The results of our covalent angled model support nonperpendicular covalent bonding sites for the DBA molecules. Because the angle and dihedral molecular degrees of freedom (eqs and 4) are involved in the expansion—much “softer” than the bond lengths (eq ), this enhances the flexibility of the framework and allows the expansion of the d001 spacing of the GOF. Although current diffraction techniques cannot verify the orientation of the DBA molecules, we suggest that the small changes to DBA–graphene bond site positions proposed in our covalent angled model are viable. The results suggest that the structure of the GOF is made up of DBA pillars randomly linked to the carbons of the graphene sheets. This implies that the GOF may not form a perfect three-dimensional network but may have some dispersion. Schaeperkoetter et al.[18] observed GOF crystallite domains sizes of ∼13.7 nm, a decrease of ∼20% from GOs, which may support this picture. A more detailed analysis of how the crystallite domain size vary with the incorporation of DBA pillars and during adsorption may help clarify this point.

In conclusion, we have simulated the adsorption of three existing and one new models of GOF-DBA and find that only the model (iv) “covalent angled GOF” is compatible with the observed swelling of the d001 spacing during supercritical adsorption of methane and xenon. This should help the development of new structured porous materials for high-performance adsorption of gases and catalytic reactions.