Influence of Structural Heterogeneity on Diffusion of CH4 and CO2 in Silicon Carbide-Derived Nanoporous Carbon.

We investigate the influence of structural heterogeneity on the transport properties of simple gases in a Hybrid Reverse Monte Carlo (HRMC) constructed model of silicon carbide-derived carbon (SiC-DC). The energy landscape of the system is determined based on free energy analysis of the atomistic model. The overall energy barriers of the system for different gases are computed along with important properties, such as Henry constant and differential enthalpy of adsorption at infinite dilution, and indicate hydrophobicity of the SiC-DC structure and its affinity for CO2 and CH4 adsorption. We also study the effect of molecular geometry, pore structure and energy heterogeneity considering different hopping scenarios for diffusion of CO2 and CH4 through ultramicropores using the Nudged Elastic Band (NEB) method. It is shown that the energy barrier of a hopping molecule is very sensitive to the shape of the pore entry. We provide evidence for the influence of structural heterogeneity on self-diffusivity of methane and carbon dioxide using molecular dynamics simulation, based on a maximum in the variation of self-diffusivity with loading. A comparison of the MD simulation results with self-diffusivities from quasi-elastic neutron scattering (QENS) measurements and, with macroscopic uptake-based low-density transport coefficients, reveals the existence of internal barriers not captured in MD simulation and QENS experiments. Nevertheless, the simulation and macroscopic uptake-based diffusion coefficients agree within a factor of 2-3, indicating that our HRMC model structure captures most of the important energy barriers affecting the transport of CH4 in the nanostructure of SiC-DC.

Introduction

Understanding transport properties of fluid molecules in microporous materials is of great importance due to intriguing features of fluidic phenomena at the nanoscale.[1−3] Many industrial and scientific applications such as gas storage and separation, petroleum refining, electrochemical energy storage, nanofluidics and materials screening require fundamental understanding of fluid properties in confined spaces. In recent years, there has been rapid development of a variety of new microporous materials such as carbide derived carbons (CDC),[4−7] carbon nanotubes,[8,9] and metal–organic frameworks (MOF),[10−13] and enhanced interest in the complexity of the fluid transport in such nanostructures.[1−3,14−17] Nevertheless, the transport properties of fluid molecules in such tight confinements are not very well understood. Established approaches such as the dusty gas model (DGM)[18] or Knudsen model have limited applications in nanoscale confinements,[19,20] while Maxwell-Stefan type models[21,22] disregard the effect of fluid inhomogeneities in mixture adsorption,[23] and more rigorous models based on the Boltzman equation rely on approximations[24] that may not be accurate in narrow nanopores.[15] On the other hand, molecular dynamics simulation can in principle provide exact results for prediction of fluid properties in confined spaces based on appropriate potential models. Nevertheless, this type of simulation can potentially suffer from inadequacy of the force field used, lack of a realistic adsorbent model and requirement of large computational resources for capturing slow diffusion processes. This is particularly true in the case of disordered microporous materials, where diffusion of gas molecules takes place at considerably slow rates. The use of realistic structures is, however, now within reach, with the development of the hybrid reverse Monte Carlo (HRMC) simulation method.[25−31]

There are relatively few diffusion studies based on realistic disordered adsorbent models in the literature. Notable among these is the work of Moore et al. who studied diffusion of argon in a disordered BPL carbon model,[26] and those reported earlier by Gubbins and co-workers[27,28,30,32] on diffusion of argon and nitrogen in structural model of saccharose based porous carbons, as well as a more recent study conducted by Nguyen and Bhatia on anomalies of water self-diffusion in the disordered structure of activated carbon.[33] A key advance achieved by using realistically constructed models of disordered materials, as opposed to idealized slit pore or crystalline models for diffusion studies of fluid molecules, is to address the significant influence of surface and structural heterogeneity on transport properties of adsorbed molecules. The effect of energy distribution arising from the micropore size distribution is already shown to be the source of system heterogeneity in activated carbon materials.[34,35] Moreover, the surface roughness of the adsorbent is known to have influence on surface diffusion of microporous solids.[36] Both of these features, that is, surface roughness and micropore size distribution, have been found to exist in the amorphous structure of silicon carbide derived carbon (SiC-DC),[25] and may be expected to be inherent to carbide derived carbons in general. It is therefore important that effects of structural heterogeneities on transport properties of fluid molecules be examined in new studies of this class of adsorbents. In very narrow pores, frictional resistance arises from both intermolecular and wall–molecule collisions, with the wall–fluid events dominating transport of adsorbate molecules.[3,23] However, due to lattice abnormalities, including surface roughness and corrugation of pore walls, existence of dead-ends and bottle necks, as well as other structural defects in disordered microporous materials, fluid molecules experience a very heterogeneous potential energy surface (PES) and a variety of different energy barriers along their diffusion paths within the solid matrix. Therefore, investigation of structural heterogeneity and the significance of the associated energy barriers is critical to the understanding of fluid transport in microporous materials.

Here, we use an experimentally validated atomistic model of disordered microporous SiC-DC developed in this laboratory,[25] to perform lengthy molecular dynamics simulations and to determine self-diffusivities of different gases over a broad range of loadings at various temperatures. The results are important to the understanding of transport properties of industrially important gases including carbon dioxide and methane in the disordered structure of SiC-DC, suggested to be promising for adsorption and gas storage applications due to its narrow pore size distribution and high pore volume.[25,37,38] We investigate the effect of structural heterogeneity on the internal resistances to the diffusion of methane and carbon dioxide in the microporous structure of SiC-DC, comparing the results from MD simulation with experimental measurements of gas diffusion using both microscopic quasi elastic neutron scattering (QENS) and macroscopic volumetric adsorption methods. We also investigate heterogeneity of the energy landscape in our carbon model, which is an important consequence of structural disorder and internal resistances, impacting diffusion rates and molecular accessibility in the system.[39,40] In addition to the use of molecular dynamics simulation and experimental methods for calculation of molecular diffusion, we analyze the free energy landscape of the system based on the computational method of Sarkisov[41,42] to provide more insights into structural heterogeneity and internal constrictions in the SiC-DC carbon, which have rate limiting effects on gas diffusion. We also employ the Nudged-Elastic Band (NEB) technique[43−45] to case study anomalies of molecular transition through ultranarrow pores.

Computational Details

We have employed an atomistic model of silicon carbide-derived carbon, developed in our laboratory,[25] based on experimental structure factor data obtained from neutron scattering using 50 nm particle size SiC-derived carbon, using the Hybrid Reverse Monte Carlo modeling technique.[31,46,47] The model (illustrated in Figure 1a–c) provides the spatial positions of 3052 carbon atoms in a 40 Å cubic simulation cell representing the disordered structure of SiC-DC. Details of the modeling technique, as well as validation procedure of the HRMC constructed model is discussed in our recent publication.[25] This model provides a solid matrix with topology and morphology consistent with that of the experimental sample and has been used for the entire simulations reported in this paper. All simulations performed for this study have made use of all-atom molecular models and force field parameters, detailed in the Supporting Information, as well as our recent work on HRMC modeling of SiC-DC structure.[25]

Figure 1

Illustration of atomistic structure of HRMC constructed model of SiC-DC.

Free Energy Landscape of Disordered System

Investigation of the energy landscape of a microporous system is crucial for understanding fluid transport properties in this class of materials. In very narrow pores, frictional resistance arises from both intermolecular and wall-molecule collisions, however due to tight confinement, diffusion of fluid molecules is predominantly affected by wall–fluid interactions.[3] Under this condition, significantly narrow pores provide highly energetic adsorption sites, which can potentially act as restrictive energy barriers during fluid diffusion. Distribution of such barriers is dictated by steric restrictions of the system associated with geometry and interconnectivity of its pore network. Hitherto, different computational methods have been employed to elucidate transport properties and diffusion behavior of fluid molecules in a variety of zeolites and MOFs, based on determination of the local or overall energy landscape of the system;[39,41,48−51] nevertheless, this type of calculation is overwhelmingly limited to noncarbonaceous materials. Here, we have employed the Helmholtz free energy analysis method of Sarkisov[41,42] to determine energy heterogeneity of the disordered SiC-DC carbon and to estimate important properties such as overall energy barriers of the system for different adsorbate molecules, as well as their Henry constant at infinite dilution. This method is a direct extension of the method developed earlier by Haldoupis et al., which performs its calculation based on determination of potential energy surface rather than Helmholtz free energy of binding.[52]

Assuming the simulation is performed at infinite dilution, atomic interactions are reduced to solid–fluid interactions. Under such conditions, one can relate the Henry constant with “the Helmholtz free energy changes associated with transferring a single molecule of adsorbate from the ideal gas system of volume Vs to a sample of porous structure of the same volume”, as described by Sarkisov.[42]In this equation, R and T are gas constant and temperature, ρs is the density of porous solid, KH is the Henry constant, A is the Helmholtz free energy of binding, and Z is the configurational integral of a single adsorbate molecule in the porous structure, while s and i.g. subscripts stand for solid and ideal gas states. The configurational integral Zs follows[42]where k is the Boltzmann’s constant and Us is the interaction potential energy between a fluid molecule and the solid structure. The integral is over all possible orientations and positions of the center of mass of the molecule within the adsorbent matrix. In practice, the simulation cell is divided into a large number of cubelets and the configurational integral is calculated for each cubelet. If the cubelet is small enough, the calculation is reduced to an integration taken solely over orientational configurations of the molecule. We have used 0.5 Å size cubelets to construct the grid network for calculation of free energy. The calculation is performed at every cubelet allowing us to obtain the free energy map of the solid structure. This approach also facilitates calculation of local properties such as local Henry constant or selectivity at different parts of the adsorbent model. To construct a plausible percolation path based on the energy map obtained from computation of free energy of bindings at every individual cubelet (A), we need to seek the threshold energy, at which one single molecule can traverse the system. This is in fact a combination of percolation path analysis[53] and free energy calculation. For a molecule traveling from one side of the simulation cell to the other side through a pattern of interlaced cubelets, interaction energies vary at every different spots. Thus, an estimation of “limiting free energy barrier” can be obtained from difference of the energy of the most attractive site along this path with the “percolating free energy threshold” explained above.[41,42]

Nudged-Elastic Band (NEB) Calculations

A central problem in diffusion studies of adsorbate molecules through confined spaces is estimation of transition rate (hopping rate) between two sides of an energy barrier within harmonic approximation of the transition state theory (hTST).[54,55] This cannot be achieved unless the activation energy barrier of the system is already known, following[56]where khTST is the transition rate constant, Esdl is the energy of the saddle point, Einit is the local potential energy minimum corresponding to the initial state, and v are the corresponding normal-mode frequencies. The determination of the potential energy maximum (Esdl) along the minimum energy path (MEP), the lowest energy pathway between two stable configurations, is important for calculation of transition rates of rare events. It also provides interesting information on magnitude of local energy barriers, as well as roughness and heterogeneity of the transition pathway.

Several computational techniques have been developed to accurately determine the activation energy barrier of a system or, in other words, to determine the highest saddle points along the MEP.[45,57] One of the most promising approaches is the so-called Nudged Elastic Band (NEB) technique, a class of “chain of states” methods in which a chain of several connected intermediate states are defined between two end points of a transition path, which are located at two local minima of the potential energy surface (PES). Spring interactions are imposed between neighboring images to generate a continuous path. These intermediate replica are then simultaneously optimized during simulation in such a way that the forces acting on the replica are minimized.[45] This algorithm employs a “force projection” or “nudging” feature to guarantee no competition between “true forces” and “spring forces”, as explained by Henkelman and Jonsson.[55] Although the NEB method is able to properly estimate position and magnitude of the saddle point between two given initial and final stable configurations, the method is likely to fail in search for the highest saddle point along MEP, considering there might be several minima and saddle points due to heterogeneity of the energy surface. As emphasized by Henkelman and co-workers,[56] in order for the NEB method to find the highest saddle point, it is essential to have an accurate estimate of the shape of the MEP. To address this issue, they have slightly modified the NEB algorithm, so that after regular NEB has converged to its MEP, the image with highest energy is selected to climb up to the top of the barrier in such a way that “the climbing image moves up the potential energy surface along the elastic band and down the potential surface perpendicular to the band”.[56]

In this work, transition path and saddle-points of two ultranarrow pore entries in the rigid HRMC constructed model of SiC-DC have been investigated using implementation of NEB in LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) package[58] based on the discussion and improvements of the algorithm by Henkelman et al.,[56] Henkelman and Jonsson,[55] and Nakano.[59] The minimization procedure has been performed using the “fire” style damped dynamics method, as described by Bitzek et al.[60] and implemented in LAMMPS. Except initial configurations of initial and final replica at the two ends of the transition path, which are explicitly defined, initial configurations of the intermediate replica are linearly interpolated between the first and the last replica. The initial configurations of the initial and final replica at the two ends of the transition path, in addition to the configurations of non-NEB fluid molecules have been obtained using Grand Canonical Monte Carlo (GCMC) simulation. The non-NEB fluid molecules have been considered as fixed molecular bodies merely to provide an appropriate background force field for the NEB molecule, while it crosses the barrier.

Molecular Dynamics Simulation

We have performed Equilibrium Molecular Dynamics (EMD) simulations for CO2 and CH4 in a periodic SiC-DC model over a wide range of loadings at 323 K, 600 and 1000 K using LAMMPS simulation package.[58] The simulations were performed in the canonical (NVT) ensemble, in which translational and rotational degrees of freedom of rigid bodies were both thermostated using the Nose–Hoover algorithm with chains, as originally described by Hoover[61] and Martyna et al.;[62,63] the rigid-body algorithm for NVT integration is explained elsewhere.[64] A Verlet time integrator was used with time step equal to 1 fs. Short-term intermolecular interactions were modeled using the 12–6 Lennard–Jones potential with a cutoff distance of 18 Å. The standard Ewald formalism was employed for electrostatic interactions with cutoff distance of 18 Å in such a way that pairwise interactions within this distance were computed directly and those outside this distance were calculated in reciprocal space. This way the cutoff distance became effectively infinite. Depending on the loading, MD simulations were run for 30 ns in average in the production phase so that displacement of the center of mass of the molecules was a multiple of the simulation cell dimension.

To calculate self-diffusivity of CO2 and CH4, mean-squared displacements (MSDs) of the center of mass of the molecules were collected in the Fickian regime, in which log–log dependence of MSD with time is linear. Self-diffusivity was then obtained using the well-known Einstein equation:where r(t) is the center of mass position vector of molecule i at time t, N is the number of molecules, and d is dimensionality of the system. From the self-diffusion coefficients obtained at various temperatures, we have also estimated the Arrhenius activation energy followingwhere D and D0 are diffusivity and temperature-independent pre-exponential factor, respectively.

Results and Discussion

Free Energy Analysis of Disordered SiC-DC

As briefly discussed in the previous section, the free energy landscape of a disordered system can be obtained at certain temperatures through calculation of the free energy of binding of a guest molecule inserted in every cubelet of a fine grid network, superimposed on a three-dimensional atomistic model of the adsorbent structure. The most useful information that can be extracted from this type of calculation is the percolating free energy threshold and limiting free energy barrier of diffusing molecules. This information is summarized in Table 1 for different adsorbate molecules in the HRMC constructed model of SiC-DC at 300 K. As presented in this table, methane and carbon dioxide are the most favorable adsorbates. These two molecules also experience the highest limiting free energy barrier, when traversing the system. Information given in Table 1 is visualized in Figure 2a for CO2 and Figure 2b for CH4, where the percolation paths of adsorbate molecules are compared with most energetically favorable pore spaces (areas with free energy of binding between −40 and −20 kJ/mol). As illustrated here, the most favorable pore spaces are widely scattered across the system and do not span a continuous percolation path. This implies that molecules adsorbed in these regions cannot easily leave the cavity to diffuse through the system unless they overcome the corresponding activation energy barriers for diffusion. Based on the information obtained from analysis of free energy landscape of the system, the position and strength of the most favorable and unfavorable pore spaces for the adsorption of a particular molecule can be identified. This information can be linked to the geometry of the pore space to provide further insight on the adsorption behavior of fluid molecules in confined spaces. Figure 3a–c visualizes repulsive and attractive pore spaces for a CO2 probe molecule across the model based on free energy of binding at 300 K. As depicted here, visual observation reveals that repulsive regions (pore spaces with very high positive energy values) consist of either extremely confined spaces inside highly narrow pores or the areas that are very close to the pore walls. On the contrary, attractive regions, which can favorably accommodate guest molecules, mainly occupy the inner part of the pore spaces.

Table 1

Characterization of Percolation Path Based on Analysis of the Free Energy of Binding at Infinite Dilution and 300 K

adsorbate	minimum free energy (kJ/mol)	percolating free energy threshold (kJ/mol)	limiting free energy barrier (kJ/mol)
argon	–27.31	–10.83	16.47
CO2	–38.83	–17.12	21.71
CH4	–33.50	–13.76	19.73
H2Oa	–20.17	–8.10	12.07

SPC/E water model has been used for water in this calculation.

Figure 2

Percolation path of (a) CO2 and (b) CH4 compared to high energy pore-spaces at 300 K. Dotted blue area represents percolation path of each molecule, while solid volumes (red or gold) refer to the high energy pore spaces.

Figure 3

Visualization of repulsive and attractive pore spaces for CO2 at 300 K, with the red volumes representing repulsive area (having free energy of binding between +1000 to +1766.8 kJ/mol) and the blue volumes representing pore spaces with negative free energy of binding (attractive interactions).

We note here that our analysis has considered a rigid structure, neglecting the effect of vibrations. In the literature there have been conflicting views on the extent of difference on transport properties due to the flexibility of the host lattice. For example, Jakobtorweihen et al.[65,66] report some reduction of the low-density transport coefficient for CH4 and C2H6 in a carbon nanotube when vibrations are considered through a boundary thermostat, with only a small reduction in the value at high densities. In contrast, Bernardi et al.[67] in their study of couette flow, find that the use of a boundary thermostat leads to a larger slip at the wall in comparison to the rigid wall. While not considered here, the effect of such vibrations will be investigated in subsequent studies.

According to the correlation between free energy of binding and the Henry constant given by eq 1, the free energy map of the system can also provide useful information on affinity of adosorbate molecules toward disordered carbon structure at infinite dilution. This information can be calculated from the Henry constant of the molecule or from interpretation of differential heat of adsorption, following[42,68,69]where ΔH is differential enthalpy of adsorption and n denotes adsorbed quantity. Values of the Henry constant for different adsorbate molecules, as well as corresponding differential heat of adsorption, are summarized in Table 2. In this table, we compare the Henry constant and differential heat of adsorption obtained from free energy analysis of the system with the Henry constant and isosteric heat of adsorption calculated from the Widom insertion method[70] and GCMC simulation at infinite dilution. Figure 4 illustrates the variation of the logarithm of the Henry constant with respect to the reciprocal temperature showing the temperature dependence of this property. Affinities of different adsorbate molecules for the adsorbent structure can be inferred from this plot. As depicted here, carbon dioxide and methane possess the highest affinity for adsorption in the disordered structure of SiC-DC, as the logarithm of their Henry constants is distinctly larger than those of argon and water at the same temperature. It is important to note that CH4 has a smaller limiting free energy barrier compared to CO2, despite its somewhat larger geometry. Moreover, adsorption of this molecule over different temperatures is less exothermic compared to carbon dioxide, indicative of lower affinity. These two properties together can potentially facilitate higher mobility of CH4 compared to CO2. This is shown to be true in the following sections; where the self-diffusivity of the two adsorbate molecules is calculated using MD simulation. The lowest affinity (the least exothermic adsorption process) according to Figure 4 belongs to water. This information in addition to the fact that minimum free energy of binding for water is much higher than other molecules across the system (given in Table 1) is an indication of hydrophobicity of the microporous structure of SiC-DC, which is consistently reported in the literature for similar carbonaceous materials.[71−74]

Table 2

Henry Constant and Differential Heat of Adsorption for Different Adsorbate Molecules

	temperature (K)	KH (mol/kg·bar), eq 1	KH (mol/kg·bar), GCMC	ΔH (kJ/mol), eq 6	ΔH (kJ/mol), eq 7
argon	100	94.917135 × 106	31.0404 × 106	–26.86	–20.28
	200	26.07	14.45	–24.52	–19.85
	300	0.64	0.45	–18.07	–15.95
	600	0.03	0.03	–14.82	–14.76
CO2	273	65.64	41.02	–38.12	–34.00
	300	16.89	10.71	–33.41	–31.66
	600	0.07	0.058	–25.93	–22.46
	1000	0.01	0.01	–21.20	–22.36
CH4	300	2.70	1.84	–26.07	–22.63
	600	0.05	0.04	–19.59	–18.41
	1000	0.01	0.01	–17.92	–19.82
H2O	300	0.17	0.15	–10.93	–11.30
	600	0.02	0.02	–10.60
	1000	0.01	0.01	–11.79

Figure 4

Variation of logarithm of the Henry constant with reciprocal temperature for different adsorbate molecules.

Local Energy Barriers and Transition Path

In this section, transition path and position of local saddle points for two nominated pore entries have been studied, solely to demonstrate how heterogeneity of pore structure, as well as molecular geometry can dictate different scenarios for diffusion path of a single molecule independent of its kinetic energy. This information can be used to calculate the local potential energy barrier at each pore entry, as opposed to the overall energy barrier of the system obtained from analysis of the free energy landscape. The former is a local potential energy barrier that an individual adsorbate molecule has to overcome, when it tries to cross a chosen pore entry, although the latter is an estimation of overall free energy barrier that a molecule experiences, while traversing the simulation cell from one side to the other. The method for calculation of local energy barrier is the Nudged Elastic Band (NEB) technique, as explained earlier in this paper. Two ultranarrow pore entries (named pore entry A and B) were selected for this calculation. Each pore entry provides a limited access to the cavity behind it through a narrow window, as illustrated in Figure 5. Pore entry A (Figure 5a) is a 13.30 Å long and 5.54 Å wide window (center to center distance). Pore entry B (Figure 5b), however, has almost a circular geometry with an average diameter of 7.0 Å. To ensure that initial configurations of the transiting molecules are closely located to the mean energy pathway at two ends of the transition path, Grand Canonical Monte Carlo (GCMC) simulation were used at similar thermodynamic conditions (∼8.2 mmol/g solid and 323 K) for both CO2 and CH4 molecules. The same temperature and pressure were applied to obtain positions of the non-NEB (background) molecules. During the NEB calculation, framework and background molecules were held fixed and only one NEB molecule at a time was allowed to move. Figure 6 depicts energy profiles of individual CO2 and CH4 molecules hopping through pore entries A and B. The potential energy in this figure is relative to the potential energy of the NEB molecule at the start of its transition path. Here, we note that for the noncircular pore entry A, the energy profile of the methane molecule is very smooth compared to that of carbon dioxide, showing only one maximum (Figure 6a). This is due to different molecular geometries of CO2 and CH4. As indicated earlier, a molecule can face multiple saddle points along its transition path, which is the case for CO2 here. As a linear molecule, CO2 is able to adjust itself according to the geometry of the pore mouth by rotating around its symmetry axis, while hopping through the window. This is shown in Figure 7, as well as a movie provided in the Supporting Information. As shown here, the O–C–O angle is nearly parallel to the pore entry cross-section at the starting and final points of the transition path (Figure 7a and c, respectively), while its orientation is perpendicular to the pore entry at the pore mouth (Figure 7b). Such adjusting rotations help the molecule to find a pathway with the smallest energy barrier, which in this case is obtained when the orientation of CO2 is perpendicular to the pore mouth, as illustrated in the snapshot in Figure 7b. In comparison, the spherical geometry of methane makes no difference among different possible orientational configurations of this molecule, thus leading a less noisy energy profile. A similar comparison for pore entry B is even more insightful (Figure 6b). In this case, not only CH4 holds a smooth diffusion energy profile, but CO2 also shows analogous behavior. Such similarity can be explained based on circular geometry of pore entry B, so that orientational movements of CO2 cause limited differences in the interaction potential energy of the molecule with its surrounding pore wall atoms.

Figure 5

Geometry of the pore entries A (a) and B (b) selected for NEB calculations.

Figure 6

Transition energy profile of CO2 and CH4 for pore entries A and B under the effect of background molecules (a, b), as well as pore entry B without any background molecule (c).

Figure 7

Transition of CO2 molecule through pore entry A at reaction coordinates (a) 0, (b) 0.73, and (c) 1. Snapshot (b) illustrates the orientation of the linear molecule around its symmetry axis at the highest saddle point in order to adjust itself with the limiting geometry of the pore mouth.

Figure 6 reveals further information on the effect of structural heterogeneity on the dynamics of fluid molecules in confined spaces. The local energy barrier for CH4 is larger than that of CO2 in pore entry A (Figure 6a), as may be expected from the larger molecular size of methane. However, our NEB simulations indicate that this may not be the only possible scenario. According to the energy profiles of CH4 and CO2 in pore entry B, CO2 may experience a larger energy barrier despite its smaller and more adjustable molecular geometry (Figure 6b). This suggests the importance of structural heterogeneity more directly. Considering CO2 adsorbs stronger (compared to CH4) due to its higher potential strength, a larger repulsive energy will be imposed on the molecule in ultranarrow pores, where atoms interact over distances smaller than their equilibrium distance, as for pore entry B. Similarly, when the NEB CO2 molecule is inside the cavity and outside the pore mouth, it will be attracted more strongly by the surrounding solid atoms, thus, having a lower level of energy compared to CH4. Consequently, a higher energy barrier for CO2 molecule is expected, since the difference between energy levels of the molecule inside the cavity and at the pore entry is larger. However, in pore entry A, where the pore mouth is larger, CO2 will face a less severe repulsive barrier.

Here, it is also important to study effect of the fluid–fluid interaction on the energy barrier. We have repeated our NEB calculations for pore entry B in the absence of background fluid molecules (Figure 6c). This way, we have been able to solely investigate interactions of the solid structure with a single NEB molecule without any interference from surrounding sources. The result signifies no dramatic change in magnitude of the energy barrier, as well as shape of the energy profile, although saddle-points are seen to be slightly shifted along the reaction coordinate axis, as illustrated in Figure 6c. According to this figure, CO2 still has to overcome a larger energy barrier compared to CH4, due to its higher potential strength and tight geometry of the pore entry as postulated above. Consequently, it is evident that in the limit of narrow pores, the magnitude of the energy barrier is largely dictated by the solid atoms at the pore mouth and the effect of fluid–fluid interactions is insignificant.

Our finding of the possibility of CO2 molecules being hindered more severely by ultramicropores, together with other evidence from analysis of the limiting free energy barrier at infinite dilution, discussed in the previous section, provides a basic explanation for the question that “why CH4 can diffuse faster compared to CO2 in microporous materials”,[75−77] especially since existence of such ultramicropores has been theoretically and experimentally confirmed in this class of materials.[25,78,79] Higher mobility of CH4 in the microporous structure of SiC-DC is demonstrated in the subsequent section of the current paper, based on the diffusivity results obtained from molecular dynamics simulation. Nevertheless, the implication of our NEB results on this issue should be considered inadequate for a definite conclusion, and need to be further investigated using methods such as Transition Path Sampling (TPS) to obtain hopping rates of the fluid molecules at the pore entry, taking into account effect of entropy on molecular diffusion.

In summary, our NEB simulations indicate that transition behavior of fluid particles, including shape of the energy profile, as well as quantity and magnitude of saddle points are very sensitive to the heterogeneity of the pore structure, arrangement of the neighboring atoms, molecular geometry, and energetic characteristic of the fluid molecule, in which heterogeneity of the pore structure plays a significant role.

Molecular Diffusion

To investigate dynamics of CO2 and CH4 in the amorphous structure of SiC-DC, the loading and temperature dependence of self-diffusivities of these gases have been studied by tracing spatial positions of molecular configurations over time trajectories from multiple MD simulations. Starting from initial configurations obtained using GCMC, MD simulations were performed at 323, 600, and 1000 K and self-diffusion coefficients were determined using eq 4. The results for the self-diffusivity of CO2 and CH4 at these temperatures are shown in Figures 8 and 9, respectively. The reported diffusion coefficients are obtained within 16% standard deviation. The error at the highest loading is 11%, while it increases to 20% at low loading.

Figure 8

Loading dependence of self-diffusivity of CO2 at three different temperatures.

Figure 9

Loading dependence of self-diffusivity of CH4 at three different temperatures.

It is known that Fickian diffusion behavior of gases can be only detected in long time runs, when molecules have escaped their initial local environment and traversed the entire lattice length.[26,80] This is associated with the “residence time”, which is the time a particle is moving around inside a “cage-like” confinement before it can leave the cage to the next repeating part of the simulation cell.[80] At very short time scales, on which intermolecular collisions are negligible, molecules are in the ballistic (free-flight) regime with a quadratic time dependence.[26,80] This regime is hard to detect in microporous materials, requiring MD simulation to be run at very small time steps. This is because of extremely limited spaces available to fluid molecules inside the micropores, so that molecular collisions occur shortly after simulation starts.

An “intermediate regime” is also detectable in nanoporous materials, where fluid molecules start colliding with each other, as well as with the pore wall in a local neighborhood.[80] At longer time steps, fluid molecules experience either a subdiffusive or a fully diffusive regime depending on the pore size and the level of confinement. A subdiffusive regime occurs when molecular diffusion is dominated by the confinement effect of a small pore space or by the traffic of molecules passing through a narrow window at high pressures.[26] A slope of 0.5–1 is usually measurable on logarithmic coordinates for a single-file subdiffusive regime, in which the molecules are unable to pass each other.[81−83] However, when molecules enter the diffusive (Fickian) regime, the slope will approach unity. Due to the highly heterogeneous structure of amorphous carbon in the SiC-DC model, lengthy MD simulations were performed to ensure the required residence time is met and the adsorbate molecules have explored the entire lattice length. Moreover, since the Fickian diffusion is at the focus of this study, the calculation of self-diffusivity has been performed in the region in which MSD is linearly related to the simulation time. Hence, MD trajectories were only collected after 1 ps of the simulation time in the production phase to avoid nonlinearity, as well as other anomalies associated with the ballistic and intermediate regimes. Nevertheless, as an exception to this convention, Figure 10 depicts the time dependence of the mean square displacement for CO2 and CH4 molecules after 0.1 ps of the simulation time at 323 K for both infinite dilution and ∼8.0 mmol/g loading. As illustrated here, a subdiffusive regime is observed at the beginning of the simulation before the slope of the MSD approaches unity. This has been constantly observed for both gases in this study, indicating the effect of confinement in the microporous structure of SiC-DC. At low loading, tight confinement of small pores dictates the subdiffusive regime; however, occurrence of this phenomenon at higher pressures suggests strong effect of pore confinement in addition to the tight molecular packing after start of the pore filling at intermediate loadings.[26]

Figure 10

Mean-squared displacement of CH4 and CO2 as a function of simulation time at 323 K, infinite dilution, and ∼8.0 mmol/g loading.

Our simulation reveals nonmonotonic behavior of molecular diffusion with a maximum in the loading dependence of the diffusion coefficient, as shown in Figures 8 and 9. This behavior is similar to diffusion of argon in disordered bituminous coal-based BPL carbon and ordered carbon replica of Faujasite Zeolite (C-FAU) studied by Gubbins and co-workers.[26] Self-diffusion coefficients of both methane and carbon dioxide are smaller at infinite dilution of lower temperatures; however, they increase at slightly higher loadings showing a maximum at intermediate densities. With further increase in the loading the self-diffusivity decreases again. Experimental evidence for diffusion of gases in active carbons and zeolite materials has shown the existence of such maxima in the concentration dependence of gas diffusivity at the level of pore filling.[84,85] PFG-NMR studies of disordered nono-porous materials have also confirmed this phenomenon.[86] Two competing mechanisms are postulated to be involved; a site-blocking entropic mechanism, as well as an energetically favored diffusion process. Slow diffusion of fluid molecules at infinite dilution is due to the existence of highly favorable adsorption sites in the amorphous structure of SiC-DC. At slightly higher concentrations most of these high affinity sites are already occupied by adsorbed molecules with increasing contribution from sites of weaker adsorption, where mobility is higher. With further increase in loading, fluid molecules fill up the bulk of the large pores, leading to significant increase in molecular collisions and steric interactions, which in turn leads to the decrease in self-diffusion coefficients.

This is very similar to what Bonilla and Bhatia have recently shown in connection with the effect of pore size distribution on diffusion in pore networks.[87] Moreover, according to Karger and Ruthven,[85] intracrystalline self-diffusivities can demonstrate different patterns of concentration dependence, ranging from constantly descending or ascending self-diffusivities to a pattern yielding a self-diffusivity variation with a maximum, which occurs in the presence of heterogeneity. Such diffusion patterns are thoroughly discussed by Keil et al. for zeolites.[86] Heterogeneity of SiC-DC structure is evident in our simulation results of the isosteric heat of adsorption, shown in Figure 11a–c. Here, wall–fluid and fluid–fluid contributions are depicted separately along with the total heat of adsorption calculated using the fluctuation formulawhere N refers to the number of molecules, U is the energy and ⟨⟩ indicates the average over the simulation run. As depicted in Figure 11a,b, the total heat of adsorption has a minimum at intermediate loadings at 323 K and a weaker one at 600 K; this indicates that isosteric heat is higher at infinite dilution due to stronger adsorption of molecules in narrow micropores, where fluid molecules are tightly packed. It is also the low loading region where molecules show smaller mobility (low diffusivity at infinite dilution). The adsorption process becomes less exothermic with increase in loading, because favorable adsorption sites are progressively occupied as loading increases. This trend steadily continuous until it reaches a minimum at intermediate loadings, and then increases again, due to stronger fluid–fluid interactions at elevated pressures, evident from the rising fluid–fluid contribution to the heat of adsorption in Figure 11. The effect of heterogeneity becomes less significant at very high temperature (1000 K), as demonstrated by the steadily increasing variation of the isosteric heat (Figure 11c) and decreasing trend of self-diffusivity with loading at this temperature (Figures 8 and 9). At 1000 K, the self-diffusivity does not show a maximum for CO2. Moreover, the maximum observed for CH4 at the intermediate loading is still lower than the self-diffusivity of this molecule at infinite dilution, and could also be due to statistical scatter of the MD. This behavior is consistent with the monotonically increasing trend of isosteric heat at this temperature, illustrated in Figure 11c. It should be noted here that isosteric heat of adsorption is always dominated by the wall–fluid contribution, which is an obvious indication of pore confinement. Fluid–fluid interactions become only slightly significant at higher pressures. Figures 8 and 9 also demonstrate that, while self-diffusion coefficients of both CO2 and CH4 increase with temperature as expected, at high loadings the effect of temperature on enhancement of self-diffusivity is much weaker in comparison to that at infinite dilution or intermediate loadings, an indication of tight molecular packing under these conditions.

Figure 11

Variation of isosteric heat of adsorption with loading for CH4 and CO2 at (a) 323, (b) 600, and (c) 1000 K.

A series of molecular dynamics simulation studies based on reconstructed models of carbonaceous materials have demonstrated similar behavior in diffusion of simple gases. Moore et al.[26] have shown a similar nonmonotonic behavior of the self-diffusion coefficient for the adsorption of argon in BPL carbon. Gubbins and co-workers observed an analogous trend for diffusion of nitrogen and argon in a Reverse Monte Carlo (RMC) constructed model of saccharose-based carbon.[27−29] Finally, Nguyen et al.[30] find similar behavior for their HRMC constructed model of saccharose-based activated carbon, based on the investigation of the self-diffusivity of argon and nitrogen.

The results from our MD simulations demonstrate anisotropic diffusion of fluid molecules within the SiC-DC pore network; so that self-diffusivities of CO2 and CH4 in the X direction are always larger than those in other directions (Figures 12 and 13). On the other hand, the self-diffusivities in the Y and Z directions are very close. Such behavior indicates anisotropic structure of amorphous SiC-DC, which may be a consequence of unidimensional propagation of the C-SiC interface during synthesis. Our finding for anisotropic diffusion of methane and carbon dioxide is in line with our previous studies on structural characterization of silicon carbide-derived carbon, demonstrating a disordered structure with percolation paths propagating anisotropically across the system.[25,37]

Figure 12

Loading dependence of self-diffusivity in different directions for CO2, at (a) 323 K, (b) 600 K and (c) 1000 K.

Figure 13

Loading dependence of self-diffusivity in different directions for CH4 at (a) 323, (b) 600, and (c) 1000 K.

Experimental Validation of Molecular Diffusion Coefficient

In this section, we report on our investigations of the diffusion of CH4 in Silicon Carbide-Derived Carbon (SiCDC), comparing microscopic molecular dynamics simulation with experimental Quasi-Elastic Neutron Scattering (QENS) and macroscopic uptake-based kinetics data at low densities.

Among a variety of experimental techniques suitable for the measurement of different ranges of the mean length of the diffusion path, the Quasi-Elastic Neutron Scattering (QENS) technique is known to be consistent with predictions of MD simulation.[88,89] While zero length column (ZLC) and uptake-rate measurements can only provide information on long-range diffusion,[90] pulsed-field gradients–nuclear magnetic resonance (PFG-NMR) offers more flexibility in measuring distances from a few hundreds of nanometers up to a few micrometers.[91,92] In contrast, the QENS technique is able to measure diffusion distances of the order of a few nanometers,[92] similar to MD simulation within current computational capabilities. The results from microscopic measurements such as QENS and molecular dynamics simulation do not account for rate limiting internal resistances and structural defects that are usually observable for the mean diffusion path at larger length scales,[90,92] thus, measuring diffusivity that is a few order of magnitudes larger than those from macroscopic or long-range methods.[90,92−98]

Volumetric uptake kinetics measurements of CH4 on the SiC-DC samples have been carried out at different temperatures using a Micromeritics ASAP 2020 adsorption analyzer, in which the transient pressure variation during small uptake steps was monitored. The corresponding uptake-time curve was interpreted using the model of diffusion in a bidisperse solid to yield a particle scale collective diffusivity, and a slow diffusivity in grain-scale ultramicropores. The details of the experimental study will be provided elsewhere. Here, the low density particle scale diffusion data is compared with the results of QENS and MD simulation, since self-and collective diffusivities of gases are identical in the limit of infinite dilution. The diffusivities used here have been obtained from uptake data at 400 mmHg using 22.8 μm diameter particles and were similar to those at lower pressures. Consequently, they are expected to represent the low density self-diffusivity of methane. In effect, the volumetric uptake-rate method measures the flux through the porous structure of adsorbent material under well-defined boundary conditions, based on the transient pressure change in the sample cell. The diffusivity is then calculated by matching the experimental flux to the solution of the Fick’s law-based diffusion equation.[90,99] Measurement of self-diffusivity using QENS experiment is, however, based on broadening in the elastic peak of the energy distribution of an incident neutron beam.[98−100] In practice, interaction of neutrons with diffusing particles gives rise to the Doppler shift, which in turn accounts for the broadening of the elastic peak.[99] In this study, the QENS experiments were carried out using the time-of-flight spectrometer IN6, at the Institut Laue-Langevin (ILL). The energy of the incident neutron beam was taken as 3.12 meV, corresponding to a wavelength of 5.12 Å. The elastic energy resolution could be fitted by a Gaussian function, whose half-width at half-maximum (HWHM) varies from 40 μeV at small wave vector transfer (Q) to 50 μeV at large Q. The SiC-DC sample was first equilibrated at 480 mbar and 300 K and the actual neutron scattering measurement was performed at constant loading (0.85 mmol/g solid) over a range of temperatures from 150–300 K. For methane, one follows the mobility of individual CH4 molecules, since the scattering is dominated by the large incoherent cross section of hydrogen.[98] Subtraction of the signal of the degassed material modifies the elastic intensity, as shown in Figure 14, because of the large small-angle scattering of the carbon. Spectra obtained at the different Q values could be fitted individually with a model consisting of isotropic diffusion, convoluted with isotropic rotation and with the instrumental resolution (Figure 14). Several spectra obtained at low Q could be fitted simultaneously using a jump diffusion model with a distribution of jump lengths. The error on the self-diffusivities is of 50%. The quality of the data was not sufficient to test anisotropic diffusion models. Figure 15 compares the temperature dependence of the diffusion coefficients for CH4 obtained from MD simulation at a loading of 1.2 mmol/g, with that from QENS at 0.85 mmol/g and macroscopic uptake-based data at low density. As seen in this figure, the QENS-based diffusion coefficients are as much as 1 order of magnitude larger than those based on the MD simulation. . This is because the length scale probed by this QENS measurement is not large enough: the lowest Q value, 0.29 Å–1, corresponds to a distance in real space of 2π/Q ≈ 22 Å. This is smaller than the size of the unit cell used in the MD calculations. However, what is remarkable here is the excellent agreement between predictions of the MD simulation with the macroscopic uptake-based data at low density. Since long-range diffusion of fluid molecules is considerably retarded by internal barriers arising from structural constrictions and disorder, macroscopic diffusivities are almost always several orders of magnitude smaller than microscopic diffusivities, usually measured by molecular dynamics simulation or QENS experiments, which probe much smaller length scales.[90,92−94,96−98,101] Interestingly, our MD results are close to the macroscopic-based diffusivities within a factor of 2–3. This kind of agreement is remarkable and indicative of the ability of our HRMC constructed model of SiC-DC to accurately capture the key internal constrictions and barriers arising from the disorder in this material, so that it can successfully reproduce both equilibrium and dynamics of gas adsorption based on the results provided here and in our previous publication.[25] Nevertheless, the higher activation energy of the macroscopic measurements, of 14.2 kJ/mol, in comparison to that from simulation, of 7.05 kJ/mol, seen in Figure 15 does indicate the presence of additional barriers not captured by the HRMC model. Further, the difference between the simulation-based activation energy and that of the uptake experiment may in part be due to the fact that experimental measurements were carried out over a narrow temperature range while our sampling in MD covers a much broader range. The self-diffusivities of CH4 and CO2 obtained from MD simulation are also compared in Figure 16, showing slightly higher activation energy of carbon dioxide (10.28 kJ/mol) compared to methane (7.05 kJ/mol).

Figure 14

Comparison between experimental and fitted QENS spectra obtained for CH4 at 300 K at different wave vector transfers: (a) 0.29, (b) 0.36, and (c) 0.41 Å–1. The negative elastic intensity is due to the subtraction of the empty SiC-DC.

Figure 15

Comparison of temperature dependence of EMD diffusion coefficient with QENS and uptake-based data, for CH4 at 400 mmHg.

Figure 16

Arrhenius plot of self-diffusivities of CH4 and CO2, obtained from EMD simulation.

Conclusion

We have investigated the heterogeneity of the energy landscape of the amorphous structure of microporous silicon carbide derived carbon using a representative atomistic HRMC model of the SiC-DC. The limiting free energy barrier, the percolating free energy threshold, Henry constant, and differential heat of adsorption of water, carbon dioxide, methane, and argon have been determined in the limit of infinite dilution using analysis of free energy landscape. Hydrophobicity of carbon structure is shown based on high minimum free energy of binding of water, as well as low heat of adsorption of this molecule compared to other adsorbates. It is shown that both methane and carbon dioxide are more favorably adsorbed compared to water. We have demonstrated how structural heterogeneity of the disordered system generates an inhomogeneous energy landscape, which in turn gives rise to the formation of high energy cluster-like adsorption sites for CH4 and CO2. The limiting free energy barrier of the system for each gas was computed from the difference between the most energetic energy sites and the percolating free energy threshold. This barrier is directly related to the existence and distribution of the above-mentioned high energy adsorption clusters, which affect mobility and diffusion of fluid molecules.

With the use of Nudge-Elastic Band (NEB) method, we have provided tangible evidence that molecular geometry, surface roughness, and structural disorder of the pore wall influence diffusion through ultranarrow pore entries. Our NEB calculations confirm that the diffusion of single molecules can be quite sensitive to structural heterogeneity and associated barriers. We find that the CO2 molecule can face a larger energy barrier compared to CH4 in ultramicropores, despite its linear and smaller molecular geometry. This is confirmed not only by our NEB calculations, but also by analysis of the free energy map of the system at infinite dilution.

We have also investigated self-diffusion of methane and carbon dioxide over a wide range of densities and temperatures using both simulation and experimental techniques. Anisotropic diffusion of CH4 and CO2 is found, and considered to be related to structural heterogeneity and pore size distribution of SiC-DC, arising from unidimensional motion of the reaction interface during chlorination of the SiC precursor. The structural heterogeneity and disorder of the SiC-DC is evident from the initial increase of self-diffusivity and decrease of the isosteric heat of adsorption with increase in loading. This observation is shown to be consistent with other investigations on heterogeneous carbonaceous materials. A comparison of our simulation results with experimental QENS and low density macroscopic uptake-based data shows remarkable agreement between MD-based and the macroscopically measured diffusivities. Such agreement indicates adequacy of our HRMC constructed model of SiC-DC in capturing the internal heterogeneity and barriers affecting transport in the structure. Nevertheless, from the difference between the activation energies obtained from MD simulation and QENS with that from macroscopic uptake, we conclude that there are some long-range internal barriers and structural constrictions, which are not captured by MD or the QENS experiment.