Computational models

All investigated zeolite structures (Supplementary Fig. ¹) were simulated under their purely silicate form and were maintained rigid during the simulation, with framework atoms fixed to their crystallographic positions. For CHA, STT, and BETAPA type zeolites the atomic positions were taken from the IZA database. For the zeolite beta, whose crystal structure is formed by an inter-growth of two crystallographic forms, designed, respectively, A and B, we have focused uniquely on the A polymorph in our simulation work. The silicalite-1 is known to exist in three distinct forms: a monoclinic one with Pnma space group and two orthorhombic forms with Pnma and P212121 space groups, designed respectively as MONO, ORTHO, and PARA. silicalite-1 is usually synthesized in the ORTHO form and passes into the low temperature MONO system after calcination. Further, this low temperature form (MONO) experiences a reversible phase transition into the ORTHO structure at about 350 K, as well as under adsorption of various molecules. Finally, upon adsorption of certain adsorbates such as nitrogen at 77 K or p-xylene the ORTHO structure moves into the PARA one. Thus, in order to reproduce in a basic manner the experimental isotherms, we have considered separately the low pressure (ORTHO) and the high pressure (PARA) systems as previously done by Snurr et al.³¹. The atomic positions for the ORTHO silicalite-1 structure were taken from the IZA database³², whereas those for the PARA silicalite-1 were taken from van Koningsveld et al.³³.

The investigated molecules are described through a “united atom” model that has been employed successfully for the investigation of their adsorption behavior in zeolites³⁴ and MOFs³⁵. In this model, each –CH_x– (with 0 ≥ x ≤ 3) and (=O) group is treated as a single interaction site. Such “united atoms” are connected by bonds maintained at fixed distances. In addition, several interaction sites bear a partial charge, contributing to the interaction energy through the Coulombic term. The nitrogen molecule is described via an explicit model: each nitrogen atom of the rigid molecule constitutes a single interaction center bearing a negative partial charge. In order to compensate the negative charges on the nitrogen atoms, there is a positive partial charge bearing no-interacting site in the middle of the nitrogen-nitrogen bond. Such three sites model allows reproducing the experimentally measured quadrupole moment of the nitrogen molecule. The intermolecular interactions between the adsorbate molecules were modeled using a sum of repulsion-dispersion potential term expressed as the Lennard–Jones interaction and the Coulombic interaction. The cross LJ terms were calculated applying the Lorentz–Berthelot combination rules. The bond distances, the partial charges and the interatomic potential parameters for all investigated molecules are summarized in Supplementary Table ³. Within the frame of the selected models, the nitrogen and p-xylene molecules are considered as rigid; therefore no intramolecular interactions are taken into account. While the intramolecular interactions for the acetone molecules are described solely via an harmonic bending term, for the n-hexane an additional dihedral torsion angle term is considered expressed by a cosine series potential. The parameters corresponding to those terms have been taken from the Transferable Potential for the Phase Equilibrium (TraPPE) forcefield, respectively, for n-hexane³⁶, p-xylene³⁷, acetone³⁸, and nitrogen³⁹ (initially fitted to reproduce the liquid/vapor coexistence curves of various fluid molecules and summarized in Supplementary Table ³).

The absolute adsorption isotherms of n-hexane, acetone, p-xylene at 298 K and nitrogen at 77 K were computed in each zeosil using the Monte Carlo simulation within the Grand Canonical ensemble implemented within the code MCCCS Towhee⁴⁰. These simulations consisted of evaluating the average number of adsorbate molecules whose chemical potential equals those of the bulk phase for given chemical potential and temperature. The chemical potential values were calculated by the test particle Widom insertion method from the NpT ensemble Monte Carlo simulation. The conventional scheme of the GCMC simulation for flexible, long chain molecules is expensive in computational time as the fraction of successful insertion moves is too low. Moreover, it does not fully explore the conformation part of the configuration space and thus does not allow achieving a proper distribution of bending and dihedral angles. The configurational-biased algorithm overcomes such shortcomings by sampling more efficiently the configuration space²³. We have applied the coupled – decoupled biased selection scheme developed by Martin and Siepmann, performing a coupled biased selection for Lennard-Jones and torsion angles selection steps, while decoupling the angle bending energy into split biased selections. The detail description of the particular algorithm employed in our simulation can be found in Martin and Siepmann^41,42. The structures of considered zeosils were treated as rigid and the periodic conditions were applied. A typical Monte Carlo run consisted of 3 × 10⁶ steps. Each step corresponded to a single MC move, including a center of mass translation, center of mass rotation, insertion of a new molecule, deletion of a randomly selected existing molecule, partial or complete regrowth of the adsorbate. The Ewald summation technique was used in the calculation of the long-range electrostatic interactions.