understanding transport 
 properties of fluid molecules in microporous 
 materials is of great importance due to intriguing features of fluidic 
 phenomena at the nanoscale .
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 ) , carbon nanotubes , and metal  organic frameworks 
 ( mof ) , and enhanced interest in the complexity of the fluid transport in 
 such nanostructures .
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 ) or 
 knudsen model have limited applications in nanoscale confinements , while maxwell - stefan type models disregard 
 the effect of fluid inhomogeneities in mixture adsorption , and more rigorous models based on the boltzman 
 equation rely on approximations that 
 may not be accurate in narrow nanopores .
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 .
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 , and those reported earlier 
 by gubbins and co - workers 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 .
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 .
moreover , the surface roughness 
 of the adsorbent is known to have influence on surface diffusion of 
 microporous solids . 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 ) , 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 .
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 , 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 .
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 .
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 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 to case study anomalies of molecular transition through 
 ultranarrow pores .
we have employed 
 an atomistic model of silicon carbide - derived 
 carbon , developed in our laboratory , 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 .
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 .
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 .
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 . 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 ; nevertheless , this type of calculation 
 is overwhelmingly limited to noncarbonaceous materials . here ,
we have 
 employed the helmholtz free energy analysis method of sarkisov 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 .
assuming the simulation is performed at infinite dilution , 
 atomic interactions are reduced to solid  fluid interactions . 
 under such conditions
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.1 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 . 

the configurational 
 integral zs follows2where 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 ( ai ) , 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 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 .
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 ) .
this can not be achieved unless 
 the activation energy barrier of the system is already known , following3where k is 
 the transition rate constant , e is 
 the energy of the saddle point , e is 
 the local potential energy minimum corresponding to the initial state , 
 and vi are the corresponding 
 normal - mode frequencies .
the determination of the potential energy 
 maximum ( e ) 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 .
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 ) . 

these intermediate replica are then simultaneously 
 optimized during simulation in such a way that the forces acting on 
 the replica are minimized .
feature to guarantee no competition between  true forces  
 and  spring forces  , as explained by henkelman and jonsson .
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 , 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  . 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 based on the discussion and improvements of the algorithm by henkelman 
 et al . ,
the 
 minimization procedure has been performed using the  fire  style damped dynamics method , as described 
 by bitzek et al . 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
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 .
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 .
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 and martyna et al .
short - term intermolecular interactions were modeled 
 using the 126 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:4where ri(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 following5where d and d0 are diffusivity and temperature - independent pre - exponential 
 factor , respectively .
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 
 can not 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 .
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 .
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 .
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 . 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 . 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 , 
 following6where 
 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 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 .
variation of logarithm of the henry constant 
 with reciprocal temperature 
 for different adsorbate molecules . 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 ) .
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 
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 .
geometry of the pore 
 entries a ( a ) and b ( b ) selected for neb calculations . 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 ) .
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 ) . 

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
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  , especially 
 since existence of such ultramicropores has been theoretically and 
 experimentally confirmed in this class of materials .
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 . 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 error at the highest 
 loading is 11% , while it increases to 20% at low loading . loading dependence 
 of self - diffusivity of co2 at three 
 different temperatures . 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 .
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 . at very short 
 time scales ,
on which intermolecular collisions are negligible , molecules 
 are in the ballistic ( free - flight ) regime with a quadratic time dependence .
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 . 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 .
a 
 slope of 0.51 is usually measurable on logarithmic coordinates 
 for a single - file subdiffusive regime , in which the molecules are 
 unable to pass each other .
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 .
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 .
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 .
pfg - nmr studies of disordered 
 nono - porous materials have also confirmed this phenomenon .
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 .
moreover , according 
 to karger and ruthven , 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 .
heterogeneity of sic - dc structure is evident in our simulation 
 results of the isosteric heat of adsorption , shown in figure 11a  c . here ,
fluid 
 contributions are depicted separately along with the total heat of 
 adsorption calculated using the fluctuation formula7where
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
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 . 

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 .
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 .
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 .
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 .
loading dependence of self - diffusivity in different directions 
 for co2 , at ( a ) 323 k , ( b ) 600 k and ( c ) 1000 k. loading dependence of self - diffusivity 
 in different directions 
 for ch4 at ( a ) 323 , ( b ) 600 , and ( c ) 1000 k. 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 . while zero length column ( zlc ) and uptake - rate measurements can 
 only provide information on long - range diffusion , 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 .
in contrast , the qens technique is able to measure diffusion distances 
 of the order of a few nanometers , 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 , thus , measuring diffusivity 
 that is a few order of magnitudes larger than those from macroscopic 
 or long - range methods .
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 .
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 . 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 . 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 150300 k. for methane , 
 one follows the mobility of individual ch4 molecules , since 
 the scattering is dominated by the large incoherent cross section 
 of hydrogen .
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 . 

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  , 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 .
interestingly , our md results 
 are close to the macroscopic - based diffusivities within a factor of 
 23 .
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 .
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 ) .
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  .
the negative 
 elastic intensity is due to the subtraction of the empty sic - dc .
comparison of temperature dependence 
 of emd diffusion coefficient 
 with qens and uptake - based data , for ch4 at 400 mmhg .
arrhenius plot of self - diffusivities 
 of ch4 and co2 , obtained from emd simulation .
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 .
 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 23 , indicating 
 that our hrmc model structure captures most of the important energy 
 barriers affecting the transport of ch4 in the nanostructure 
 of sic - dc . 