Sorption, Structure and Dynamics of CO2 and Ethane in Silicalite at High Pressure: A Combined Monte Carlo and Molecular Dynamics Simulation Study

Silicalite is an important nanoporous material that finds applications in several industries, including gas separation and catalysis. While the sorption, structure, and dynamics of several molecules confined in the pores of silicalite have been reported, most of these studies have been restricted to low pressures. Here we report a comparative study of sorption, structure, and dynamics of CO2 and ethane in silicalite at high pressures (up to 100 bar) using a combination of Monte Carlo (MC) and molecular dynamics (MD) simulations. The behavior of the two fluids is studied in terms of the simulated sorption isotherms, the positional and orientational distribution of sorbed molecules in silicalite, and their translational diffusion, vibrational spectra, and rotational motion. Both CO2 and ethane are found to exhibit orientational ordering in silicalite pores; however, at high pressures, while CO2 prefers to reside in the channel intersections, ethane molecules reside mostly in the sinusoidal channels. While CO2 exhibits a higher self-diffusion coefficient than ethane at low pressures, at high pressures, it becomes slower than ethane. Both CO2 and ethane exhibit rotational motion at two time scales. At both time scales, the rotational motion of ethane is faster. The differences observed here in the behavior of CO2 and ethane in silicalite pores can be seen as a consequence of an interplay of the kinetic diameter of the two molecules and the quadrupole moment of CO2.


Introduction
The behavior of molecules is strongly determined by their environment, and is different if they are surrounded by molecules of the same species or of a different species. The latter scenario occurs, for example, when a fluid is confined in a porous media, or is adsorbed on the surface of a substrate. The specific interaction between the fluid molecules and the confining substrate alters the behavior of the fluid. As most of the interesting chemical reactions occur at interfaces between two species, the behavior of molecules under confinement or at an interface in general is an interesting topic of research [1][2][3][4][5]. An important component that determines the interaction between molecules in general is the distribution of the electrostatic charge. In some molecules such as water, an asymmetrical distribution of charge gives rise to a permanent dipole moment. The distribution of charge in a molecule plays an important role in determining the dynamical behavior. For example, it has been found that for molecules of similar shapes and sizes, both the translational and the rotational motions are severely affected by an asymmetry in the charge distribution [6,7]. The motion of acetonitrile, acetaldehyde, and acetone, all with a non-zero dipole moment was found to be severely restricted, whereas that of propane with no dipole moment was found to be less restricted. In the studies mentioned above, the motion of these molecules confined in all silica analogues of ZSM-5-silicalite was studied. Silicalite imposes strict geometrical confinements on the sorbed molecules, due to its small channel-like pores, which are, on average, 0.55 nm in diameter. With the kinetic diameters of the sorbed molecules being of a size that is comparable to the pore diameter in these studies, an interplay between the effects of the dipole moment and the geometrical confinement on these molecules resulted in severely restricted motion for all molecules except for propane. CO 2 is a linear molecule with no dipole moment, but with a finite quadrupole moment. This quadrupole moment of CO 2 gives rise to interesting phenomena. For example, it has been suggested that due to a stronger interaction with quadrupolar CO 2 , this molecule replaces hydrocarbons from silica pore walls. This has been shown to give rise to an enhancement in the diffusivity of confined hydrocarbon by the addition of CO 2 in several studies [8][9][10]. Juxtaposing the quadrupolar nature of CO 2 with the effects of dipole moments on the motion of molecules in silicalite reported in [6,7], a question can be raised whether the quadrupolar nature of CO 2 can also result in a severely restricted motion of CO 2 in silicalite. To address this, a detailed study of CO 2 dynamics in silicalite is needed. The sorption and mobility of gases such as CO 2 and ethane inside ZSM-5 has been studied by several authors [11][12][13][14][15]. More recently, mixtures containing CO 2 in silicalite have also been studied [16][17][18]. Moisture-driven CO 2 sorption has also been studied. For example, Shi et al. have reported the effects of moisture on the sorption of CO 2 in nano-structured sorbents [19,20]. However, most of these studies were limited to lower loadings corresponding to pressures of a few tens of bar. Studies at higher pressures are lacking. In particular, there are very few studies addressing molecular motion under confinement at supercritical densities. The effects of high densities encountered in nanoporous engineered and natural materials are thus underrepresented in the literature.
In molecular level computer simulations, the force field used is an important parameter that can determine the course of the simulation. These force-fields are usually obtained by fitting some experimental data. For small hydrocarbons and CO 2 , the TraPPE (transferable potentials for phase equilibria) force field developed by Martin and Siepmann [21] is commonly used. Alkanes are generally modeled in united atom (UA) formalism in this force field, and an ethane molecule is modeled as a dumbbell-type diatomic molecule with two CH 3 pseudo-atoms. An ethane molecule modeled in the TraPPE-UA force field is thus similar to the linear CO 2 molecule, although with a kinetic diameter of 0.44 nm [22], the ethane molecule is slightly bigger than the CO 2 molecule with a kinetic diameter of 0.33 nm [23] (the kinetic diameter is a measure of the sphere of influence of a molecule, and it represents the size of a molecule in reference to intermolecular collision events). A comparative study of CO 2 and ethane can reveal the role of the quadrupole moment in the sorption and dynamics of these two molecules inside silicalite. Here, we report grand canonical Monte Carlo (GCMC) and molecular dynamics (MD) studies on the sorption, structure, and dynamics of ethane and CO 2 in silicalite at higher loadings corresponding to a partial pressure of up to 100 bar. To achieve supercritical densities for both CO 2 as well as ethane, this study was carried out at a temperature of 308 K. At this temperature, the bulk state of CO 2 and ethane becomes supercritical at 46 bar and 77 bar, respectively [24]. Figure 1a shows the sorption isotherms of ethane and CO 2 in silicalite at 308 K. At pressures below atmospheric pressure (1 bar), ethane is adsorbed in silicalite more favorably. However, at higher pressures, the sorption of CO 2 is more favorable. The supercritical pressures of bulk ethane and CO 2 are marked in the figure by using vertical dashed lines of black and red colors, respectively. The density of a fluid adsorbed in a nanoporous material at a given temperature and pressure condition is generally higher than the bulk density at those conditions. It is therefore customary to express the sorption isotherms in terms of the excess density (n ex ) of the sorbed fluid over the corresponding bulk density [25]. The relation between the sorbed density (n ads ) and the excess density can be expressed as: where V p is the volume of the pores accessible to the fluid, M is its molar mass, and ρ b is the bulk density of the fluid at the corresponding temperature and pressure conditions. The quantity n ads is obtained directly from the GCMC simulation. First et al. [26] used an automated approach to characterize several zeolitic frameworks, and obtained pore volumes. For our calculations of the excess sorption, we used their estimate of the accessible pore volume in silicalite. The values of the bulk density of several fluids including CO 2 and ethane are available on the NIST web-book [24]. Using these quantities, the excess sorption isotherms of ethane and CO 2 obtained at 308 K are shown in Figure 1b. Both ethane and CO 2 show a region of enrichment in the pores, up to a pressure of 20 bar, followed by depletion. In what follows, we shall discuss the results from MD simulations for ethane and CO 2 in silicalite, carried out on loadings corresponding to three subcritical and one supercritical pressure (0.1, 1, 20 and 100 bar). These pressures are also marked in Figure 1a.

Sorption Isotherms of CO 2 and Ethane in Silicalite
where Vp is the volume of the pores accessible to the fluid, M is its molar mass, and ρb is the bulk density of the fluid at the corresponding temperature and pressure conditions. The quantity nads is obtained directly from the GCMC simulation. First et al. [26] used an automated approach to characterize several zeolitic frameworks, and obtained pore volumes. For our calculations of the excess sorption, we used their estimate of the accessible pore volume in silicalite. The values of the bulk density of several fluids including CO2 and ethane are available on the NIST web-book [24]. Using these quantities, the excess sorption isotherms of ethane and CO2 obtained at 308 K are shown in Figure 1b. Both ethane and CO2 show a region of enrichment in the pores, up to a pressure of 20 bar, followed by depletion. In what follows, we shall discuss the results from MD simulations for ethane and CO2 in silicalite, carried out on loadings corresponding to three subcritical and one supercritical pressure (0.1, 1, 20 and 100 bar). These pressures are also marked in Figure 1a.

Distribution of the Molecules in Different Pores
Silicalite has a complex pore structure, with pores in the form of straight and sinusoidal channels, with their intersections exhibiting somewhat larger dimensions then the channels themselves. While the straight channels are oriented along the Cartesian Y direction, the sinusoidal channels lie in the X-Z plane. Due to the slightly different shapes of these pores, we explored which of these three pore types affords hosting of the guest molecules more readily. To probe this, we obtained the distribution of the molecules amongst these three pore types, in all simulations. These distributions are shown in Figure 2 for both ethane and CO2, with black and red symbols, respectively. At lower pressures, the difference in the number of ethane molecules is similar in all channels, within one standard deviation. However, at higher pressures, ethane molecules show a clear preference for residing in the sinusoidal channels. This could be because a molecule of large enough kinetic diameter might find it difficult to exit out of sinusoidal channels, due to its tortuosity. On the other hand, CO2 has a smaller kinetic diameter, resulting in an increased ability to move out of these channels, and therefore it shows no preference to reside in them, rather preferring intersections.

Distribution of the Molecules in Different Pores
Silicalite has a complex pore structure, with pores in the form of straight and sinusoidal channels, with their intersections exhibiting somewhat larger dimensions then the channels themselves. While the straight channels are oriented along the Cartesian Y direction, the sinusoidal channels lie in the X-Z plane. Due to the slightly different shapes of these pores, we explored which of these three pore types affords hosting of the guest molecules more readily. To probe this, we obtained the distribution of the molecules amongst these three pore types, in all simulations. These distributions are shown in Figure 2 for both ethane and CO 2 , with black and red symbols, respectively. At lower pressures, the difference in the number of ethane molecules is similar in all channels, within one standard deviation. However, at higher pressures, ethane molecules show a clear preference for residing in the sinusoidal channels. This could be because a molecule of large enough kinetic diameter might find it difficult to exit out of sinusoidal channels, due to its tortuosity. On the other hand, CO 2 has a smaller kinetic diameter, resulting in an increased ability to move out of these channels, and therefore it shows no preference to reside in them, rather preferring intersections.

Orientational Distribution of Molecules in Silicalite and Anisotropy
In a matrix like silicalite, that puts severe geometrical constraints on the sorbed molecules, both positional as well orientational ordering of molecules can be expected. Indeed, orientational ordering has been observed for several molecules that have been sorbed in silicalite [6,15,27]. Orientational ordering is observed overall for both ethane and CO2 in silicalite (See Supplementary Materials, Figure S1). Due to the tortuous shape of the sinusoidal channels a wide distribution of orientations can be expected for confined molecules. Likewise, at channel intersections, the geometrical restriction is slightly relaxed, due to a larger pore volume and ellipsoidal shape; thus, orientational ordering of sorbed molecules is expected to be less constrained, resulting in relatively more disorder. In straight channels that are parallel to the Cartesian Y axis, the ordering is relatively stronger, and the distribution is easier to interpret. Therefore, in Figure 3a, we show the orientational distribution function (ODF) of ethane and CO2 at 100 bar in the straight channels of silicalite, with respect to the three Cartesian directions. In absence of a preferred orientation, the distribution is expected to be isotropic, and the corresponding ODF is also included in Figure 3a for comparison. Any deviation from the isotropic ODF signifies a preferred orientation. Clearly, both ethane and CO2 show deviations from isotropic behavior. However, the preferred directions in the two cases are different. While ethane prefers to lie mostly along the channel axis (with the Y-direction showing stronger than isotropic distribution at lower (closer to 0°) and higher (closer to 180°) angles), CO2 prefers to lie perpendicular to the channel axis (with the Y-direction ODF peaking at 90°, and being weaker at lower and higher angles, compared to isotropic ODF).
The extent of the deviation in the orientational distribution with respect to the isotropic case can be used to quantify the anisotropy in the orientational structure of the sorbed molecules. Following previous studies [6,15], we define the degree of anisotropy (ϕ) in the orientational structure as: where gi(θ) is the orientational distribution function, the subscript i denotes the case of isotropic distribution (iso) or the silicalite (s), and the sum is taken over N = 3150 different values of θ over which the ODF were calculated. The degree of anisotropy that is defined in this manner is shown in Figure 3b for the three Cartesian directions for ethane and CO2 in straight channels of silicalite at all pressures. For ethane, the variation of ϕ with pressure shows a consistent trend in all directions; namely, ϕ decreases with an increase in pressure. For CO2, however, the trend is different in different

Orientational Distribution of Molecules in Silicalite and Anisotropy
In a matrix like silicalite, that puts severe geometrical constraints on the sorbed molecules, both positional as well orientational ordering of molecules can be expected. Indeed, orientational ordering has been observed for several molecules that have been sorbed in silicalite [6,15,27]. Orientational ordering is observed overall for both ethane and CO 2 in silicalite (See Supplementary Materials, Figure S1). Due to the tortuous shape of the sinusoidal channels a wide distribution of orientations can be expected for confined molecules. Likewise, at channel intersections, the geometrical restriction is slightly relaxed, due to a larger pore volume and ellipsoidal shape; thus, orientational ordering of sorbed molecules is expected to be less constrained, resulting in relatively more disorder. In straight channels that are parallel to the Cartesian Y axis, the ordering is relatively stronger, and the distribution is easier to interpret. Therefore, in Figure 3a, we show the orientational distribution function (ODF) of ethane and CO 2 at 100 bar in the straight channels of silicalite, with respect to the three Cartesian directions. In absence of a preferred orientation, the distribution is expected to be isotropic, and the corresponding ODF is also included in Figure 3a for comparison. Any deviation from the isotropic ODF signifies a preferred orientation. Clearly, both ethane and CO 2 show deviations from isotropic behavior. However, the preferred directions in the two cases are different. While ethane prefers to lie mostly along the channel axis (with the Y-direction showing stronger than isotropic distribution at lower (closer to 0 • ) and higher (closer to 180 • ) angles), CO 2 prefers to lie perpendicular to the channel axis (with the Y-direction ODF peaking at 90 • , and being weaker at lower and higher angles, compared to isotropic ODF). directions, with ϕ decreasing at higher pressures monotonously in the Z and X directions, while the variation in the Y direction is non-monotonous.

Translational Motion and its Anisotropy
Translational motion of ethane and CO2 in silicalite is studied by using the mean squared displacement (MSD) of the sorbed molecules. In an anisotropic matrix such as silicalite, MSD in different directions can be expected to be different. In Figure 4a, we show the MSD of ethane and CO2 in the three Cartesian directions X, Y, and Z, along with the overall MSD at 100 bar. At this pressure, ethane can be seen to exhibit consistently higher MSD than CO2 in all directions. This is in spite of the smaller kinetic diameter of CO2, as compared to ethane. A significant degree of anisotropy can be observed in the MSD, with the Y direction showing the highest MSD for both ethane and CO2, whereas the Z direction shows the smallest MSD. This is because of the straight channels which because of their geometry, facilitate the motion of the sorbed molecules that lie along the Cartesian Y directuion, whereas the sinusoidal channels with their tortuous shape inhibiting motion, lie in the X-Z plane. MSD is a measure of the self-diffusive motion of molecules and at long times, it can be used to obtain the self-diffusion coefficient (Ds) using the Einstein relation: where nd is the number of degrees of freedom, and t is time. For consistency, we used the slope of MSD vs. t curve between t = 400 ps and t = 600 ps, to obtain the self-diffusion coefficient in all cases.
In Figure 4b, we show the self-diffusion coefficient of ethane and CO2 (solid black and red symbols respectively) as a function of pressure. The variation of Ds with pressure is similar for both ethane and CO2. Ds decreases at higher pressure, as can be expected, due to the crowding of molecules at higher pressures. At lower pressures, ethane exhibits slower diffusion compared to CO2, whereas at higher pressures, ethane becomes relatively faster. This could be because of the relatively lower amounts of ethane sorbed in silicalite at higher pressures. In Figure 4b, we also show the ratio of selfdiffusivity in the directions Y and Z (Ds Y /Ds Z ; open symbols) as a function of pressure. This ratio is a measure of the anisotropy in the translational diffusion. The anisotropy in translational diffusion is found to increase at higher pressures for both ethane and CO2. Further, the anisotropy at higher pressures is greater for ethane, compared to CO2. This could be a consequence of the slightly larger The extent of the deviation in the orientational distribution with respect to the isotropic case can be used to quantify the anisotropy in the orientational structure of the sorbed molecules. Following previous studies [6,15], we define the degree of anisotropy (φ) in the orientational structure as: where g i (θ) is the orientational distribution function, the subscript i denotes the case of isotropic distribution (iso) or the silicalite (s), and the sum is taken over N = 3150 different values of θ over which the ODF were calculated. The degree of anisotropy that is defined in this manner is shown in Figure 3b for the three Cartesian directions for ethane and CO 2 in straight channels of silicalite at all pressures. For ethane, the variation of φ with pressure shows a consistent trend in all directions; namely, φ decreases with an increase in pressure. For CO 2 , however, the trend is different in different directions, with φ decreasing at higher pressures monotonously in the Z and X directions, while the variation in the Y direction is non-monotonous.

Translational Motion and its Anisotropy
Translational motion of ethane and CO 2 in silicalite is studied by using the mean squared displacement (MSD) of the sorbed molecules. In an anisotropic matrix such as silicalite, MSD in different directions can be expected to be different. In Figure 4a, we show the MSD of ethane and CO 2 in the three Cartesian directions X, Y, and Z, along with the overall MSD at 100 bar. At this pressure, ethane can be seen to exhibit consistently higher MSD than CO 2 in all directions. This is in spite of the smaller kinetic diameter of CO 2 , as compared to ethane. A significant degree of anisotropy can be observed in the MSD, with the Y direction showing the highest MSD for both ethane and CO 2 , whereas the Z direction shows the smallest MSD. This is because of the straight channels which because of their geometry, facilitate the motion of the sorbed molecules that lie along the Cartesian Y directuion, whereas the sinusoidal channels with their tortuous shape inhibiting motion, lie in the X-Z plane. MSD is a measure of the self-diffusive motion of molecules and at long times, it can be used to obtain the self-diffusion coefficient (D s ) using the Einstein relation: 6 of 13 where n d is the number of degrees of freedom, and t is time. For consistency, we used the slope of MSD vs. t curve between t = 400 ps and t = 600 ps, to obtain the self-diffusion coefficient in all cases.
In Figure 4b, we show the self-diffusion coefficient of ethane and CO 2 (solid black and red symbols respectively) as a function of pressure. The variation of D s with pressure is similar for both ethane and CO 2 . D s decreases at higher pressure, as can be expected, due to the crowding of molecules at higher pressures. At lower pressures, ethane exhibits slower diffusion compared to CO 2 , whereas at higher pressures, ethane becomes relatively faster. This could be because of the relatively lower amounts of ethane sorbed in silicalite at higher pressures. In Figure 4b, we also show the ratio of self-diffusivity in the directions Y and Z (D s Y /D s Z ; open symbols) as a function of pressure. This ratio is a measure of the anisotropy in the translational diffusion. The anisotropy in translational diffusion is found to increase at higher pressures for both ethane and CO 2 . Further, the anisotropy at higher pressures is greater for ethane, compared to CO 2 . This could be a consequence of the slightly larger kinetic diameter of ethane, which makes it more difficult for ethane molecules to diffuse through sinusoidal channels (X-Z plane), compared to the straight channels.  Another important quantity that can be used to understand the translational motion of molecules is the velocity autocorrelation function (VACF, Cv(t)).

Cv(t) = <v(t + t0).v(t0)>
where v(t) is the velocity of an atom at time t. The angular brackets denote an average over the number of atoms, as well as the time origins t0. The self-diffusion coefficient can also be obtained by integrating the Cv(t) vs. t curve [28]. In addition, Cv(t) can be used to calculate the vibrational spectrum (I(ω)) of the system, using the relation [29]: We point out that since both ethane and CO2 are modeled as rigid molecules here, the vibrational spectra obtained using the above equation will be limited to the intermolecular vibrational modes alone. In Figure 5, we show the vibrational spectra calculated from the VACF for ethane (top panel) and CO2 (bottom panel) at different pressures. While the spectra for ethane are bimodal at low pressures, the lower energy peak merges with the broad peak at higher pressures. CO2 on the other hand exhibits a single broad peak at all pressures. Further, a shift to higher energies can be observed in the spectra for both ethane and CO2 at higher pressures. A similar shift to higher energies has been observed in the vibrational spectrum of propane confined in MCM-41-S in an inelastic neutron scattering experimental study [29]. For comparison, the vibrational spectra of bulk ethane [30] and CO2 [31] obtained from inelastic neutron scattering experiments are also included. These spectra are taken from the ISIS INS database [32]. Although the spectra calculated from a simulation of rigid molecules is not expected to capture all the details of the experimental spectra, some salient features Another important quantity that can be used to understand the translational motion of molecules is the velocity autocorrelation function (VACF, C v (t)).
where v(t) is the velocity of an atom at time t. The angular brackets denote an average over the number of atoms, as well as the time origins t 0 . The self-diffusion coefficient can also be obtained by integrating the C v (t) vs. t curve [28]. In addition, C v (t) can be used to calculate the vibrational spectrum (I(ω)) of the system, using the relation [29]: We point out that since both ethane and CO 2 are modeled as rigid molecules here, the vibrational spectra obtained using the above equation will be limited to the intermolecular vibrational modes alone. In Figure 5, we show the vibrational spectra calculated from the VACF for ethane (top panel) and CO 2 (bottom panel) at different pressures. While the spectra for ethane are bimodal at low pressures, the lower energy peak merges with the broad peak at higher pressures. CO 2 on the other hand exhibits a single broad peak at all pressures. Further, a shift to higher energies can be observed in the spectra for both ethane and CO 2 at higher pressures. A similar shift to higher energies has been observed in the vibrational spectrum of propane confined in MCM-41-S in an inelastic neutron scattering experimental study [29]. For comparison, the vibrational spectra of bulk ethane [30] and CO 2 [31] obtained from inelastic neutron scattering experiments are also included. These spectra are taken from the ISIS INS database [32]. Although the spectra calculated from a simulation of rigid molecules is not expected to capture all the details of the experimental spectra, some salient features are reproduced. For example, the intermolecular vibrational modes of CO 2 confined in silicalite peak at lower energies compared to ethane, a result that is qualitatively consistent with the experimental data on bulk ethane and CO 2 . Further, the calculated spectra of both ethane and CO 2 look more disordered, as compared to the bulk spectra. A similar confinement-induced disorder was observed in the vibrational spectra of propane in MCM-41-S [29] obtained from INS experiments, as well as MD simulations.  [30] and CO2 [31], obtained in inelastic neutron scattering experiments and available at the ISIS INS database [32].

Rotational Motion
The preferred orientations seen in the ODF in Figure 2 for both ethane and CO2 suggest that the rotational motion of the sorbed molecules is restricted. The rotational motion of rigid molecules in general can be studied by following the orientation of a unit vector (u) that is attached to the molecule as a function of time [33,34]. In particular, one can follow the evolution of the first-order rotational correlation function (RCF, CR(t)), defined as: In Figure 6a we show the RCF of ethane (black curve) and CO2 (red curve) at 100 bar. Compared to ethane, the RCF of CO2 takes longer to decay to zero, indicating a slower rotational motion of CO2 compared to ethane. Further, the RCF of both ethane and CO2 exhibit distinct behaviors at times that are shorter than a picosecond, and longer times. This facilitates dividing the RCF in two regimes of sub-picosecond times, and longer times. The boundary of these two regimes is marked clearly in case of ethane in the form of a kink. As mentioned in previous publications, this is a signature of vibrational motion [15,27]. No kink can be observed in the case of CO2 suggesting that the rotational motion of CO2 is different from the libration type of motion exhibited by ethane and other molecules [6,15,27]. CO2 has a kinetic diameter that is considerably smaller than the channel diameter of silicalite. Further, the interaction of the molecule with the pore walls is not strong enough to restrict it from rotating more freely inside the silicalite channels, tracing out the entire orientational space available.

Rotational Motion
The preferred orientations seen in the ODF in Figure 2 for both ethane and CO 2 suggest that the rotational motion of the sorbed molecules is restricted. The rotational motion of rigid molecules in general can be studied by following the orientation of a unit vector (u) that is attached to the molecule as a function of time [33,34]. In particular, one can follow the evolution of the first-order rotational correlation function (RCF, C R (t)), defined as: In Figure 6a we show the RCF of ethane (black curve) and CO 2 (red curve) at 100 bar. Compared to ethane, the RCF of CO 2 takes longer to decay to zero, indicating a slower rotational motion of CO 2 compared to ethane. Further, the RCF of both ethane and CO 2 exhibit distinct behaviors at times that are shorter than a picosecond, and longer times. This facilitates dividing the RCF in two regimes of sub-picosecond times, and longer times. The boundary of these two regimes is marked clearly in case of ethane in the form of a kink. As mentioned in previous publications, this is a signature of vibrational motion [15,27]. No kink can be observed in the case of CO 2 suggesting that the rotational motion of CO 2 is different from the libration type of motion exhibited by ethane and other molecules [6,15,27]. CO 2 has a kinetic diameter that is considerably smaller than the channel diameter of silicalite. Further, the interaction of the molecule with the pore walls is not strong enough to restrict it from rotating more freely inside the silicalite channels, tracing out the entire orientational space available.
Molecules 2018, 23, x FOR PEER REVIEW 8 of 12 when pressure is increased to 1 bar, and for higher pressures, it slows down again. For ethane, the short time rotational motion is independent of pressure, while the longer time rotation becomes faster consistently with an increase in pressure.
(a) (b) Figure 6. Rotational motion of ethane (black) and CO2 (red) in silicalite: (a) Rotational correlation function (RCF) at 100 bar. (b) Time scales of rotational motion obtained by fitting the two regimes of RCF with an exponential decay. In the first regime, τ1 was obtained by fitting the RCF between 0 and 0.2 ps. The Y-axis labels for τ1 can be read on the right in blue, while those for τ2 appear on the left.

Discussion
The sorption isotherms presented in Figure 1 are in very good agreement with those reported by Sun et al. [13] for ethane in silicalite at 308 K up to 3.43 bar, and CO2 in silicalite up to 1.37 bar. Ethane shows a stronger sorption at lower pressures, while CO2 is adsorbed more readily at higher pressures. In a study on CO2 sorption in metal organic frameworks at 298 K up to 35 bar, Walton et al. [35] found that a higher amount of CO2 is adsorbed in the metal organic framework, if electrostatic interactions are accounted for in the simulation. Conversely, lower amounts of CO2 were found to be adsorbed if the electrostatic interactions were switched off. The higher sorption of CO2 is thus due to the electrostatic interactions between the quadrupole moment of CO2 and the pore surface.
Both ethane and CO2 exhibit type I isotherms up to 100 bar. It is interesting to note that all three types of pores in silicalite intersections, straight channels, and sinusoidal channels are accessible to both ethane, as well as CO2 at all pressures, as shown in Figure 2. Even at lower pressures not investigated with MD, GCMC data show that the molecules start to populate all three types of pores (see Supplementary Materials, Figure S2). This is in contrast to the case of branched alkanes, e.g., isobutane, that populate only the intersections at lower pressures, resulting in a type IV isotherm, as observed by Vlugt et al. [36]. This is because of the relatively smaller kinetic diameter of ethane and CO2, compared to iso-butane. It is natural that CO2 shows a clear preference for intersections at higher pressure. This is because the intersections offer a marginally larger ellipsoidal free space, compared to the straight and sinusoidal channels. However, ethane molecules at higher pressures prefer to stay in sinusoidal channels, although at lower pressures, this preference is not very clear. We note that in our previous study of ethane in silicalite, the highest loading investigated (eight molecules per unit cell) corresponded to a pressure of slightly less than 1 bar. At that pressure, the error bars in Figure  2 make no conclusion about the preference for pore type difficult.
The orientational order exhibited by both ethane and CO2 (Figure 3) is consistent with that exhibited by several other molecules confined in silicalite [6,15,27]. This behavior results from a strict geometrical restriction imposed by the narrow silicalite pores. Further, the preferred orientation of CO2 in straight channels, perpendicular to the channel direction, is a consequence of the partial charges on the extremities of the molecule. The electrostatic interactions between the extremities of CO2 molecule and the pore wall make the molecule lie perpendicular to the channel, whereas in the absence of such partial charges, ethane molecules prefer a parallel orientation. This has important consequences for the motion of these molecules. Being linear molecules, it can be expected that the  Figure 6b for ethane (black) and CO 2 (red) at different pressures. For consistency, all RCFs were fit between t = 0 and t = 0.2 ps; and t = 2 ps and t = 200 ps to obtain τ 1 and τ 2 , respectively. The rotational motion of CO 2 is slower than that of ethane in both regimes. Although the short time rotational motion of CO 2 shows a monotonous decrease at higher pressures, the longer time rotation becomes faster initially, when pressure is increased to 1 bar, and for higher pressures, it slows down again. For ethane, the short time rotational motion is independent of pressure, while the longer time rotation becomes faster consistently with an increase in pressure.

Discussion
The sorption isotherms presented in Figure 1 are in very good agreement with those reported by Sun et al. [13] for ethane in silicalite at 308 K up to 3.43 bar, and CO 2 in silicalite up to 1.37 bar. Ethane shows a stronger sorption at lower pressures, while CO 2 is adsorbed more readily at higher pressures. In a study on CO 2 sorption in metal organic frameworks at 298 K up to 35 bar, Walton et al. [35] found that a higher amount of CO 2 is adsorbed in the metal organic framework, if electrostatic interactions are accounted for in the simulation. Conversely, lower amounts of CO 2 were found to be adsorbed if the electrostatic interactions were switched off. The higher sorption of CO 2 is thus due to the electrostatic interactions between the quadrupole moment of CO 2 and the pore surface.
Both ethane and CO 2 exhibit type I isotherms up to 100 bar. It is interesting to note that all three types of pores in silicalite intersections, straight channels, and sinusoidal channels are accessible to both ethane, as well as CO 2 at all pressures, as shown in Figure 2. Even at lower pressures not investigated with MD, GCMC data show that the molecules start to populate all three types of pores (see Supplementary Materials, Figure S2). This is in contrast to the case of branched alkanes, e.g., iso-butane, that populate only the intersections at lower pressures, resulting in a type IV isotherm, as observed by Vlugt et al. [36]. This is because of the relatively smaller kinetic diameter of ethane and CO 2 , compared to iso-butane. It is natural that CO 2 shows a clear preference for intersections at higher pressure. This is because the intersections offer a marginally larger ellipsoidal free space, compared to the straight and sinusoidal channels. However, ethane molecules at higher pressures prefer to stay in sinusoidal channels, although at lower pressures, this preference is not very clear. We note that in our previous study of ethane in silicalite, the highest loading investigated (eight molecules per unit cell) corresponded to a pressure of slightly less than 1 bar. At that pressure, the error bars in Figure 2 make no conclusion about the preference for pore type difficult.
The orientational order exhibited by both ethane and CO 2 (Figure 3) is consistent with that exhibited by several other molecules confined in silicalite [6,15,27]. This behavior results from a strict geometrical restriction imposed by the narrow silicalite pores. Further, the preferred orientation of CO 2 in straight channels, perpendicular to the channel direction, is a consequence of the partial charges on the extremities of the molecule. The electrostatic interactions between the extremities of CO 2 molecule and the pore wall make the molecule lie perpendicular to the channel, whereas in the absence of such partial charges, ethane molecules prefer a parallel orientation. This has important consequences for the motion of these molecules. Being linear molecules, it can be expected that the motions of both these molecules along their molecular axes will be easier. Indeed, for ethane, it was found that the MSD parallel to the molecular axis is considerably higher than that which is perpendicular to the axis (see Supplementary Materials, Figure S3). However, for CO 2 , no considerable difference was observed for MSD that was parallel or perpendicular to the molecular axis. Although, in a free environment, the motion of a CO 2 molecule along the molecular axis would be easier, in case of silicalite, no free space is available in that direction, as the molecule lies perpendicular to the channel axis.
At lower pressures, CO 2 diffuses faster than ethane. However, at those pressures, the amount of CO 2 adsorbed in silicalite is lower than that of ethane. Thus, crowding effects on the translation diffusion for ethane are stronger. However, at higher pressures, the relative amount of CO 2 is larger than ethane, and hence the diffusivity of CO 2 is slower. The variation of self-diffusivity can be viewed as the interplay of the kinetic diameters and crowding effects. The rotational motion of CO 2 in silicalite is slower than ethane, in spite of a slightly smaller kinetic diameter. While for ethane, the rotational motion gets consistently faster at higher pressures, for CO 2 the variation is non-monotonous. In a previous work on ethane rotation in silicalite [15], we observed a correlation between the variation of rotational motion and anisotropy in the ODF. We note the same correlation here in the case of ethane, wherein as the anisotropy decreases at higher pressure, the rotational motion becomes faster. A similar correlation can also be seen between the rotational motion of CO 2 and anisotropy in ODF in the Y direction.
Although the behavior of both ethane and CO 2 in silicalite at pressures corresponding to bulk supercritical pressures is quantitatively different from that at lower pressures, no major qualitative changes were observed at this pressure (100 bar) compared to the lower pressures corresponding to bulk sub-critical pressures.
In the absence of electrostatic interactions, CO 2 , which is a smaller molecule, can be expected to experience relatively less restriction to motion by the confining media. However, due to the quadrupole moment of CO 2 that appears as partial charges on the C and O atoms in the TraPPE model used here, CO 2 is not only adsorbed in silicalite to a greater extent, its motion is also suppressed more, compared to ethane. On the other hand, the suppression in the mobility of CO 2 is not as drastic as what is observed in the case of molecules with dipole moments that barely move inside silicalite [6,7]. Thus, although not as strongly as the dipole moment of other molecule, the quadrupole moment of CO 2 does affect its sorption, structure, and dynamics in silicalite.

Materials and Methods
In this work we used Grand Canonical Monte Carlo (GCMC) simulations to obtain the sorption isotherms of ethane and CO 2 in silicalite, and Molecular Dynamics (MD) simulations to study their structural and dynamical properties. While GCMC simulations were carried out using DL_Monte [37], MD simulations were carried out using DL_Poly [38]. Both these packages have been developed by CCP Forge, and they are available free of charge. DL_Monte has been shown to have a computational performance that is comparable to several other Monte Carlo codes [39]. Besides, having commonalities with DL_Poly, it was a natural choice to complement the MD simulations carried out using DL_Poly.
A model ZSM-5 pore network was built by making use of the co-ordinates provided by Koningsveld et al. [40]. To start with, two molecules each of CO 2 or ethane were loaded in the straight channels of ZSM-5. GCMC simulations were then carried out on the simulation cell thus obtained. During the simulation, the sorbed molecules, i.e., ethane or CO 2 , could be inserted, deleted, translated, or rotated, while all silicalite atoms were kept rigid. All simulations were carried out using a fixed gas partial pressure at 308 K. This is done by using the following selection procedures for the insertion (P i ) or deletion (P d ) of a molecule.
where V is the simulation cell volume, P is the partial pressure of the gas, U represents the potential energy, N is the number of molecules and β = k B T.
To start with, GCMC simulations were carried out at the lowest pressure of 0.05 bar. After this, successive simulations were carried out by sequentially increasing the pressure. The final configurations of each simulation were used as the initial configurations for the next simulation. At lower pressures, 2 million Monte Carlo steps were enough to obtain statistically meaningful configurations. Of these, the first 500,000 steps were discarded to account for equilibrium. At higher pressures, 5 million steps were used, out of which, the first 150,000 were equilibration steps. Coordinates were sampled every 10,000 steps. During the simulation, the translation, rotation, insertion, and deletion of ethane or CO 2 were undertaken with a probability of 0.25, 0.25, 0.5, and 0.5 respectively. A potential cut-off of 1.4 nm, as suggested for TraPPE-UA force field was used. The number of molecules of ethane or CO 2 in the simulation cell averaged over the entire production run, discarding the equilibration steps, were used to obtain the adsorption isotherms as a function of pressure.
Having obtained the adsorption isotherms, four pressures of interest (0.1, 1, 20, and 100 bar) were selected to further investigate the structures and dynamics of ethane and CO 2 in silicalite, using MD simulations. The outputs from the MC simulations were used as the inputs for the MD simulations in all cases. MD simulations were carried out in NVT ensemble, using a Nose-Hoover thermostat. A total of eight simulations, four simulations corresponding to the selected pressure each for ethane and CO 2 in silicalite, were carried out by using a simulation time step of 1 fs. In each simulation, the system was allowed to equilibrate at 308 K for 0.5 ns. After this, trajectories were recorded at each 0.02 ps for a duration of 1.5 ns. Various properties of interest were then calculated from these trajectories. The achievement of equilibrium was confirmed when the energy and temperature values showed no clear monotonous increase/decrease in time, and the variations were less than 5%.

Conclusions
We studied the sorption, structure, and dynamics of ethane and CO 2 in silicalite-two similar molecules without and with a quadrupole moment. The quadrupole moment of CO 2 plays an important role in all of the properties investigated here. Combining the information obtained here with some previous studies, we can conclude that the effects of electrostatic interactions on the properties of a confined fluid in a nanoporous medium are strongest for molecules with a net dipole moment, and become progressively milder for molecules with a quadrupole moment and no polarity.