Diffusion of Aromatics in Silicalite-1: Experimental and Theoretical

Subscriber access provided by Northern Illinois University

Article

Diffusion of Aromatics in Silicalite-1: Experimental and Theoretical Evidence of Entropic Barriers Panagiotis D. Kolokathis, Gyorgy Kali, Herve Jobic, and Doros N. Theodorou J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b05462 • Publication Date (Web): 18 Aug 2016 Downloaded from http://pubs.acs.org on August 21, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Diffusion of Aromatics in Silicalite-1: Experimental and Theoretical Evidence of Entropic Barriers Panagiotis D. Kolokathis*†, György Káli‡,¶, Hervé Jobic§, and Doros N. Theodorou*†

†

School of Chemical Engineering, National Technical University of Athens, Zografou Campus, GR-15780 Athens, Greece. ‡ ¶

Institut Laue-Langevin, 6 rue Jules Horowitz, 38000 Grenoble, France

Wigner Research Centre for Physics, Konkoly Thege Miklós út 29-33, H-1121 Budapest, XII

§

Institut de Recherches sur la Catalyse et l’ Environnement de Lyon, CNRS, 2 av. Albert Einstein, 69626 Villeurbanne, France,

*

E-mail: [email protected], [email protected]. Telephone: +30 210 772 3157

Abstract A specific computational methodology based on transition state theory (Kolokathis, P.D. et al., Mol. Sim., 2014, 40, 80–100) is evolved and applied for the calculation of the self-diffusion coefficients of para-xylene and benzene in silicalite-1 at infinite dilution. In addition, we study the orientational distributions of phenyl rings and methyl stems of para-xylene and benzene sorbed in the zeolite and check for the existence of entropic barriers to translational motion. A new reduction method for the states appearing in the free energy profiles is presented and used for the calculation of transition rate constants for elementary jumps. Quasi-Elastic neutron scattering measurements are also conducted and compared with the simulation results. A major conclusion from both experiments and simulations is that para-xylene diffuses roughly 100 times faster than benzene when sorbed at low occupancy in silicalite. Benzene encounters strong entropic barriers to translational motion at channel intersections, where it can adopt a variety of orientations. The corresponding barriers for para-xylene are much lower, reflecting the inability of its major axis to reorient within channel intersections.

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 42

I.

INTRODUCTION

The study of diffusion of aromatics in silicalite-1 is a topic of major importance for the petroleum industry.1 Silicalite-1 channels (straight and sinusoidal) have approximately 5.5Å diameter and, as a result, silicalite-1 can be used as a molecular sieve for aromatics with different kinetic diameter.2,3 A number of studies have been conducted and a large number of different experimental techniques4-17 have been used to measure diffusion coefficients of aromatics in silicalite-1. However, a controversy among these experimental techniques still exists concerning the diffusion of para-xylene (p-xylene) at infinite dilution in this crystal. Some studies4-7 show that the diffusivity of p-xylene is of the same order as that of benzene8-13 (actually slightly faster), while others14-17 support that p-xylene’s diffusivity is 2 orders of magnitude higher than that of benzene. These experiments have been performed for sorbates at infinite dilution in the Ortho18 phase of silicalite-1. Among various measurement methods, Quasi-Elastic Neutron Scattering (QENS) offers the advantage that it measures displacements over sub-microscopic dimensions and is therefore free of limitations due to internal and external mass transport resistances caused by defects, which are known to otherwise exist in this type of materials.19,20 Molecular simulations can also be used as a tool to explain these differences and to predict diffusion coefficients, as has been done in the past.21-27 For molecules experiencing a close fit within the zeolite channels, such as aromatics in silicalite-1, diffusion is too slow to be predicted reliably by molecular dynamics (MD) simulations and strategies based on infrequent event analysis of elementary jumps between sorption sites within the channels are necessary. In recent studies, a new method, MESoRReD, for the calculation of self-diffusion in crystals with orthorhombic symmetry was invented.28 This method leads to a very fast calculation of the diffusion coefficient once the rate constants of the elementary jumps of the sorbate have been calculated and the network of transitions inside the lattice is also known. Then, from the lowest nonzero eigenvalue of a matrix containing these rate constants for a single unit cell of the material, one can derive the self-diffusivity tensor. A new method for the calculation of rate constants of elementary jumps from MD simulations was also invented.27 This method relies on a modification of umbrella sampling;29 it will be referred to as Curvilinear Umbrella Sampling (CUS) in the following. To our knowledge, this was the first time dimensionality reduction of a free energy profile has been accomplished for curvilinear paths using a specific projection method and rules. In the work of Kolokathis et al.27 it has been explained that the projection definition does not allow using classical umbrella sampling along curvilinear paths. This problem does not arise for straight paths. Straight paths and classical umbrella sampling have been recently used by Verploegh et al.30 and Camp et al.31 to study the diffusion of hydrocarbon molecules in zeolitic imidazolate frameworks and light gas diffusion in a porous organic cage crystal 3 (CC3), respectively. Curvilinear paths are impossible to avoid in some crystals, however, as in the case of silicalite-1 because of its sinusoidal channels. In addition, curved paths help us describe well the dimensionality reduction of a multi-dimensional dividing surface to a single point along a one-dimensional curvilinear coordinate without losing any information. Such loss of information is common when we use straight paths. The calculation of transmission coefficients is one way to compensate for this loss, without being able to eliminate it, as can be understood from the 2-dimensional free energy ACS Paragon Plus Environment

Page 3 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


profiles of ref 27. (see also Figure S1 in the Supporting Information to the present paper). Using curvilinear paths derived from the 3-dimensional free energy landscape, one can see that transmission coefficient values depend only on the strength of coupling to the thermostat process of the crystal and, as a result, are expected to be near unity for sharply peaked free energy profiles. For this reason we have decided to apply the CUS algorithm described in ref 27 and not classical umbrella sampling in the present work. Both classical umbrella sampling and CUS27 are biased dynamics methods, while the bluemoon ensemble32 (another method used in conjunction with transition state theory) is a constraint dynamics one.33 The bluemoon ensemble cannot be used for curved paths with the definition of projection we suggest, because it is not possible to calculate the Jacobian determinant it needs. The reader is reminded that reaction coordinate-free methods, such as Transition Path Sampling (TPS),34 could have been applied to the problem. TPS offers the advantage, in comparison with conventional umbrella sampling and transition state theory-based techniques, that the dividing surface and the transmission coefficient are not needed to calculate the rate constants with accuracy.35 However, TPS is much more involved and computationally demanding than transition state theory-based methods. TPS has been used successfully in the past to study the diffusion of linear hydrocarbons in silica zeolite LTA.35 It has also been implemented for isobutane in silicalite-1; in that case, however, it yielded large error bars (uncertainty of two orders of magnitude) for the final diffusivity, due to the error introduced in matching histograms for the free energy calculation.36 According to Van Erp et al.,37 a drawback of the rate constants as calculated by TPS is that they can be highly oscillatory because of recrossings and will reach a plateau only after a relatively long time. The path length in TPS must exceed the typical timescale of these oscillations, introducing a need for long runs.37 The CUS approach adopted here takes advantage of the known crystal geometry of silicalite-1, relies on conventional MD simulations, provides useful information concerning the 3-D free energy profile associated with translational motion, readily allows the separation of energetic from entropic barriers to motion along the channels, and involves modest computational cost. This is why it was preferred to rigorous, more general, but more sophisticated and computationally demanding methods, such as TPS. An alternative to CUS would be to consider a rectilinear transition path along each channel segment. Umbrella sampling conducted using rectilinear transition paths and our CUS method should give exactly the same transition rates, if all dynamical correction factors were correctly calculated. However, CUS has the advantage of yielding a continuous 1-D free energy profile, from which all rate constants can be accurately obtained with near-unity transmission coefficients. A more elaborate discussion of this issue is given in Section I of the Supporting Information, where the pitfalls associated with projecting the full configuration space onto onedimensional transition paths are discussed. From Section I of the Supporting Information one can see that the method of ref 27 is the most suitable method for our problem.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 42

II. MOLECULAR SIMULATION DETAILS a) Sensitivity of free energy profiles to i) electrostatics calculation and ii) flexibility of the sorbate In our previous work27 we used a static electrostatic map to describe the electrostatic interactions between silicalite-1 and benzene in order to enhance the speed of the calculations. This means that partial charges were considered fixed at the crystallographic positions corresponding to mechanical equilibrium of the pure crystal,38 while the atoms of the crystal were able to move around these positions as dictated by a flexible model force field. The work of Clark and Snurr39 showed that it is necessary for the atoms of the crystal to be in their crystallographic positions, otherwise the predicted Henry’s law constant will deviate significantly from experimental values. Having this in mind, the model of Vlugt and Schenk40 was used to describe the flexibility of silicalite-1 and to keep its atoms near their crystallographic positions. According to the observations of Clark and Snurr,39 the Vlugt and Schenk model is expected to predict the Henry constant correctly. To be absolutely sure about this, a Monte Carlo method utilizing a large number of sorbate insertions (compare ref 41) should be used in conjunction with a large number of instantaneous configurations of the flexible silicalite-1. This calculation is expected to entail high computational cost in order to converge with low error bars, and was not undertaken here. On the other hand, isosteric heats of sorption have been calculated with the model of Vlugt and Schenk and found to be in very good agreement with experiment (see Section IIe below). In the present work we have decided to check the deviation in free energy profiles when we use a static electrostatic map versus explicit calculation of all electrostatics from moving partial charges on the zeolite. The latter explicit calculation is not trivial if we want to keep the atoms around their crystallographic positions, and special computational provisions are needed. The details of this explicit calculation are outlined in section II of the Supporting Information. From Figure 1 one can see that molecular simulation using a static electrostatic map is 8 times faster in 1 CPU than a simulation with moving charges on the zeolite which uses the explicit P3M method,42-44 even though the former was conducted on a larger number of atoms (3×3×3 silicalite unit cells instead of 2×2×3 of the P3M calculation). The static electrostatic map calculation is still faster, even though we use 8 CPUs in the P3M calculation. At the same time the free energy profiles are practically the same, as one can see from Figure 2. A large number of MD simulations are needed to obtain fragments of the free energy profile. It is preferable to conduct many serial (1 CPU) simultaneous simulations with static charges than parallel ones with moving charges, as Figure 1 shows that, in this way, we can achieve significant savings in computer time.


Page 5 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Computation time required a) for the system benzene-27 silicalite unit cells using a static map for the electrostatics and a time-step of δt=0.5fs (red color); b) for the system benzene-12 unit cells (blue color) with explicit summation of electrostatic interactions from moving partial charges on the zeolite and a time-step of δt=1fs. All calculations conducted on Dell PowerEdge 1950 QUAD Core Servers with Intel XEON E5430 processors at 2.66 GHz and 16GB memory. It must be noted that the subroutine which projects the center of mass of the sorbate molecule has been written for one only processor. However, the calculation of electrostatic interactions is by far the limiting step in all computations. In addition to the simulations of the rigid benzene in flexible silicalite-1 with a) static electrostatic map and b) explicit calculation of electrostatic interactions with moving charges on the atoms, we decided to study the influence of the benzene flexibility on the free energy profiles. To do this, we used the Condensed-phase Optimized Potentials for Atomistic Simulation Studies (COMPASS) model to describe intra-benzene interactions.45 Because of the high computational requirements of the P3M calculation (see Figure 1), 1-dimensional free energy fragments were calculated for the range 11.6 to 15.1 Å of the ξ coordinate of the straight channel. This range was chosen because of the large slope of free energy profiles observed inside it, which make it the most sensitive to flexibility issues. The free energy profiles calculated from the previous simulations are illustrated in Figure 2. From these results one can see that neither the use of a static electrostatic map nor the flexibility of benzene influence the 1-dimensional free energy profile. In addition, combining Figure 2 with Figure 1, one can see that the use of a static electrostatic map and a rigid benzene is the most profitable method for simulating the sorption and diffusion of benzene inside silicalite. Similar behavior is expected for the simulation of sorbates in crystals at such loadings of sorbate that a phase transition of the crystal does not occur. The method of using a static electrostatic map along with Lennard-Jones interactions which follow the oscillations of the atoms of the crystal can be very useful in the simulation of systems with very large unit cells, such as MIL100,41 especially when one examines the diffusion of molecules with similar kinetic diameter as the MIL100 window. Clark and Snurr39 have shown that electrostatics cause large differences in the average potential energy, even for small deviations of the silicalite-1 atoms from their crystallographic positions. In our calculations, however, free energy profiles obtained using ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 42

explicit electrostatics and free energy profiles based on a static electrostatic map turn out to be almost the same (Figure 2). Consequently, a force field in which atoms, as Lennard-Jones centers, oscillate around their correct crystallographic positions but Coulombic interactions associated with partial charges are handled through a static electrostatic map can be considered satisfactory.

Figure 2. a) Local 1-dimensional free energy profile of benzene in silicalite-1 in the straight channel for a temperature of 555Κ using FF1(analytical calculation of electrostatic interactions with a flexible zeolite framework and a rigid benzene), FF2(analytical calculation of electrostatic interactions with a flexible zeolite framework and a flexible benzene), FF3(static electrostatic map and rigid benzene). b) free energy fragments of the area between the walls (where there is no influence of the bias potential of the wall), from which the 1-dimensional local free energy profile has been derived by the stitching process.  is the reaction coordinate along the straight channel. The labels S and I indicate the sorption sites (free energy minimum regions) corresponding to the benzene center of mass residing in the straight channel and in the channel intersection region, respectively


Page 7 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


b) Choice of force field for the description of intra-silicalite-1 interactions The force field of Vlugt and Schenk40 is a modified edition of the force field of Demontis et al.46 Actually, the Vlugt and Schenk40 force field uses the crystallographic equilibrium distances for all the bonds, while the Demontis et al.46 force field uses specific mean equilibrium values. As mentioned above, Clark and Snurr39 found it necessary to keep the zeolite atoms in their crystallographically correct positions38 in order to predict the Henry’s law constant and the isosteric heat in close agreement with experimental values. The Demontis et al.46 force field describes the IR spectra of zeolites well; the same is expected of the Vlugt and Schenk40 force field, as they use the same spring constants. However, the Vlugt and Schenk40 force field fails to keep the volume of silicalite-1 (the lattice actually collapses) when used in MD simulations under constant pressure (Parrinnello-Rahman isothermal-isostress algorithm with characteristic time τP=1000fs)47,48 and constant temperature Nosé-Hoover algorithm with characteristic thermostat relaxation time τΤ=100fs)48,49 simulation, and so it can be used only in constant volume simulations. The above characteristic times are defined by Melchionna et al.48 in the isothermal-isostress algorithm they present. From now on, we use the term DHFF-O (Demontis Harmonic Force Field Original) for the Demontis et al. force field and the term DHFF for the Vlugt and Schenk one. A different force field for zeolites has been proposed by Smirnov and Bougeard50 (we use the term SGVFF-O (Simplified General Valence Force Field Original) for this one). They introduced harmonic angle potentials in addition to a harmonic bond potential for the Si-O bond. They also used specific mean equilibrium values for all the bond lengths and angles. We need to modify this force field by using the equilibrium crystallographic bonds and angles38 for the reasons mentioned before, if we want to study the diffusion of aromatics in silicalite-1. We call this modified force field SGVFF. From our calculations, by conducting constant pressure and temperature simulations with the parameters we mentioned, we calculated a bulk modulus 20.63GPa for the SGVFF-O and 40.46GPa for the SGVFF. We remind the reader that the experimental value is 18.2 GPa.51 Elastic constants, such as the bulk modulus, provide a measure of dimensional changes of the crystal lattice resulting from the imposition of external mechanical stress. A force field providing reasonable estimates of the elastic constants would likely capture internal distortions of the lattice brought about by the sorption of bulky, tight-fitting molecules such as the aromatics. These distortions are of paramount importance in shaping the free energy barriers at the bottlenecks to translational motion encountered by these molecules. Too high (low) values of the bulk modulus would indicate that the crystal model is too stiff (soft) and therefore that the free energy barriers to translational motion might be over(under)estimated. The ability of the force field to capture thermal expansion of the lattice is also important, since we are interested in understanding and predicting diffusion over a range of temperatures. The unit cell volume of pure Silicalite-1 as a function of temperature, as predicted by the force fields we use, is shown in Figure 4. Figures S14-S16 of the Supporting Information display the individual dimensions of the unit cell along the a, b, and c crystallographic axes. Clearly, the negative signs of the linear thermal expansivities along the three crystallographic axes, and of the volumetric thermal expansivity, are correctly predicted.52,53 The magnitude of the volumetric



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 42

expansivity is reasonable compared to experiment, even though the predicted unit cell volumes exhibit a small systematic deviation from experimental values. In the literature one can find other force fields for silicates.52,54-60 However, these force fields incorporate explicit electrostatic and Lennard-Jones interactions among crystal atoms instead of the harmonic potentials invoked in DHFF and SGVFF. We have decided to use only the DHFF and the SGVFF force fields in our simulations of benzene and p-xylene in silicalite-1 because of their simplicity and relative success in reproducing pure silicalite-1 properties. If we use a force field cast only in terms of bond stretching and bond angle bending harmonic potentials with equilibrium bonds and bond angles equal to the ones determined from crystallography, the minimum of the total potential energy is guaranteed to be attained at a configuration coinciding with the crystallographic one. This means that the most probable configuration of the crystal in the course of a simulation will coincide with the set of crystallographic coordinates. If electrostatic interactions were introduced in addition to the bond and bond angle harmonic potentials, as happens in most force fields in the literature, then the potential energy minimum would shift, and the equilibrium positions would not be identical with the crystallographic ones. We wish to stay as close as possible to the correct crystallographic positions because, as mentioned above, Clark and Snurr39 have shown that this is important for the correct prediction of sorption properties. The comparative study of DHFF and SGVFF helps us examine how much the silicalite-1 force field choice influences our final results (free energy profiles, rate constants, diffusion coefficients). All the force fields parameters used in this work are provided in Tables S1-S7 in the Supporting Information. Here, it is important to note that harmonic potentials (bonds, angles, dihedrals, impropers) are good to be used when we are sure that the crystal does not undergo a phase transition for the loading, temperature and pressure conditions we study. For example, Fröhlich et al. recently showed that a harmonic dihedral does not let the CAU10H crystal go through a phase transition, whereas experiments show that it does.61


Page 9 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Pressure as a function of volume of silicalite-1, after 1ns of NPT MD simulation a) for SGVFF-O force field and temperature 300Κ (blue) b) for SGVFF-O force field and temperarure 465Κ (green) c) for SGVFF force field and temperature 300Κ (red) d) light orange point shows the experimental volume of Olson et al.38 e) purple star points show the volumetric behavior of silicalite-1 as synthesized using a bulk material dissolution technique51 (fluoride route) f) dark orange points show the properties of silicalite-1 which has been synthesized using NaOH, tetrapropylammonium bromide and tetraethyl orthosilicate as described in ref 51 (alkaline route).

Figure 4. Volume of silicalite-1 unit cell as a function of temperature, a) after 10ns of NPT MD simulation using SGVFF-O force field b) after 10ns of NPT MD simulation using SGVFF force field c) experimental values by Krokidas et al.52 d) experimental values from Bhange and Ramaswamy.53 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 42

c) Investigation of dynamics of benzene in silicalite-1 using SGVFF In this investigation we have used the parameters of ref 27, just replacing the DHFF with the SGVFF force field in describing intra-silicalite interactions. A time step δt=0.5fs was employed in the Verlet algorithm tracking the motion of silicalite-1 atoms in conjunction with velocity rescaling on silicalite’s atoms to keep the temperature constant.62,63 In this way, silicalite-1 works as a thermal bath for the rigid sorbate. In every molecular dynamics simulation the center of mass velocity and angular velocity distribution were calculated in addition to the velocities of silicon and oxygen atoms of silicalite-1. They all follow the expected Maxwell-Boltzmann distribution  mv 2 I  2   m  2 kBTx  I xx  2 xxkBTx    (vx )    where ρ is the probability density, vx is the ,  (x )   e e   2 k BT  2 k BT      velocity in the x direction, m is the mass, Ixx is the moment of inertia tensor for rotation around the x axis, T is the temperature and kB is the Boltzmann constant, and similarly for y and z. We picked simple velocity rescaling for thermostating our crystal, rather than a more sophisticated extended ensemble method such as Nosé-Hoover,48,49,63 because for a large number of atoms (as in our case) velocity rescaling describes well the canonical ensemble (see Figures S10-S13 in Supporting Information).64 For tracking the motion of the rigid sorbate the algorithm proposed by 63,65 Fincham was used. The electrostatics was described by the same static electrostatic map in conjuction with cubic Hermite interpolation as in ref 27. The wall potential parameters used in accumulating the center of mass distribution within different regions of the zeolite in the context of umbrella sampling were the same as in ref 27. The paths on which the benzene center of mass was projected were also the same. All these parameters are listed in more detail in Sections IV and XI of the Supporting Information. Next we conducted a large number of MD simulations, in order to get all the 1-dimensional free energy profile fragments we need and accumulate the final 1-dimensional free energy profiles for each path and for three different temperatures. These profiles and their fragments are shown in Figures 5 and 6. The fact that the free energy returns to the same value, as dictated by the periodicity of the crystal, after stitching the different fragments together, provides evidence for the consistency of our calculations.


Page 11 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. a) Free energy profiles for benzene, in RT units, as functions of the reaction coordinate at temperatures of i) 300K (blue) ii) 465Κ (red) iii) 555Κ (green) along the straight silicalite-1 channel for a crystal described by the SGVFF force field and b) free energy fragments of the area between the walls (where there is no influence of the bias potential of the wall), from which the 1-dimensional local free energy profile has been derived by the stitching process. S and I mark sorption sites (free energy minimum regions) in the straight channels and intersections.

Figure 6. a) Free energy profiles for benzene, in RT units, as a function of the reaction coordinate at temperatures of i) 300K (blue) ii) 465Κ (red) iii) 555Κ (green) along the sinusoidal silicalite-1 channel for a crystal described by SGVFF force field and b) free energy fragments of the area between the walls (where there is no influence of the bias potential of the wall), from which the 1-dimensional local free energy profile has been derived by the stitching process. Z and I mark sorption sites (free energy minimum regions) in the sinusoidal channels and intersections. ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 42

d) Investigation of the dynamics of p-xylene in silicalite-1 using DHFF and SGVFF To simulate p-xylene in silicalite-1 we have used the same algorithms (Verlet – velocity rescaling thermostat on the zeolite atoms, Fincham for the rigid penetrant),62,63,65 timesteps, wall parameters, and reaction coordinate paths as in Part IV of the Supporting Information for benzene. The parameters of the rigid p-xylene (charges, bond lengths, Lennard-Jones parameters) are described analytically in Section IV of the Supporting Information; they have been used successfully in the past to describe Henry’s constants, isosteric heats, and sorption isotherms of p-xylene in rigid silicalite-1.66, 67 The methyl group of p-xylene is represented in a united atom approximation. In contrast to the benzene case, the 3-dimensional (3-D) free energy profiles of p-xylene in silicalite-1 have not been studied in the past. We have chosen to use the paths proven to be satisfactory for the study of benzene in silicalite-1 and, by moving the repulsive walls along them, to acquire the 3-D free energy fragments for p-xylene as well. As pointed out in our past work,27 these fragments are not accurate enough for computing transition rate constants reliably because of their limited sample size and correspondingly low resolution; their main purpose is to show us the locations of the minima (states) and the dividing surfaces. The 3-D free energy profiles of p-xylene, compared with those calculated in the past in ref 27 for benzene, are illustrated in Figure 7a. Clearly, p-xylene is more restricted in the intersection than benzene. From these profiles we can see that the straight line used in ref 27 is also satisfactory as a reaction coordinate (i.e., describes the dividing surface well) for projecting the center of mass of p-xylene and ultimately determining its 1-D free energy profile. Because we have used the same path (x=10 Å, y=y, z=0.0), we can compare the 1-D free energy profiles of benzene and pxylene directly, as we can see in Figure 7b. There is a large difference between the 1-D free energy profiles in the region of the intersection. To explain this, we use Figure 7c. In this Figure, the orientational free energy profiles of the phenyl rings of benzene and p-xylene are compared. We can see that in the intersection the benzene ring can orient easily in a variety of directions (cyan color on the surfaces of the little spheres), while it actually prefers one orientation (red color) for each position of its center of mass. P-xylene, on the other hand, behaves differently. The normal vector to its phenyl ring exhibits a far less diffuse orientational distribution and a strong preference for a specific orientation (sharp red points on the surfaces of the little spheres) at each position of its center of mass. It seems to be very difficult for p-xylene to turn to other directions and its long axis seems restricted to be roughly parallel to the axis of the straight channel of the zeolite. This means that benzene in the intersection adopts a lot more microstates than does p-xylene and, as a result, its entropy in the intersection region is higher and its free energy displays a deeper minimum in the intersection region. Thus, the difference in the 1-D free energy profiles between benzene (blue) and p-xylene (red) in Figure 7b is attributable to entropic effects associated with the ability of the molecules to adopt different orientations in the intersection region. We will elaborate on this topic later. Similar conclusions can be reached if we compare the methyl orientational free energy profile of Figure 7c (and, for a better view, Figure 8a ), with the benzene hydrogen orientational free energy profile, Figure 8b. The benzene hydrogen orientations have been measured by Quasi Elastic Neutron Scaterring68 and Nuclear Magnetic Resonance in the past69 and are in perfect agreement with the orientational free energy profiles of Figure 8b. It is also important to add that the ring orientational free energy profiles of ACS Paragon Plus Environment

Page 13 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


benzene and p-xylene in the interior of the straight channel (far from the intersection) are very similar and this is the reason for the small difference between their 1-D free energy profiles within the channel. Thus, we suspect the presence of a significantly higher barrier to translational motion along the straight channel in the case of benzene relative to p-xylene. This barrier is of entropic origin and arises in the intersection region.

Figure 7. a) 2-dimensional cross-section of the 3-D free energy profile of p-xylene (left) and benzene (right) at plane x=10Å b) 1-dimensional free energy profiles of p-xylene (red) and benzene (blue) and the 1-dimensional free energy fragments from which they have been constructed c) Orientational free energy profiles of i) normal vector to benzene ring ii) normal vector to p-xylene ring and iii) carbon-methyl bond (methyl stem) of p-xylene. The orientational free energy profiles have been determined from the natural logarithm of the probability density distribution of the vector shown in red in the small molecular diagrams and are depicted as contour plots on the surface of a sphere for various values ξ of the reaction coordinate along the straight channel. In these contour plots, red color signifies high probability for the orientation considered. All calculations here have been conducted with the DHFF force field for a temperature of 300K.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 42

Figure 8. a) Orientation of one methyl of p-xylene along the straight channel of silicalite b) Orientation of one hydrogen of benzene along the same channel for a temperature of 300K. Colorbar is in RT units. All calculations have been performed with the DHFF force field.

We have performed simulations for the temperatures of 405K and 485K in addition to 300K, because, as we will present later, Quasi Elastic Neutron Scattering measurements have also been conducted at those temperatures. The 1-D free energy profiles for these higher temperatures are illustrated in Figure 9. Until now we have examined the free energy profiles in the straight channels. Next, we focus on the behavior of p-xylene in the sinusoidal channel. Similarly to the straight channel, we chose the sinusoidal channel path that has proved to be satisfactory for the study of benzene in silicalite-1. By moving the repulsive walls we have acquired the 3-D free energy fragments. We have stitched these 3-D free energy fragments together and the final 3-D free energy profile is illustrated in Figure 10.


Page 15 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9. i) Free energy profiles of p-xylene, in RT units, as a function of the reaction coordinate for 300K (blue), 405Κ (red), and 485Κ(green) along the straight silicalite-1 channel for a crystal described by the DHFF force field (a) and the SGVFF force field (b) and ii) free energy fragments of the area between the walls (where there is no influence of the bias potential of the walls), from which the 1-dimensional local free energy profile has been derived by the stitching process The path we used for the sinusoidal channel in our previous work for benzene27 and in Figure 6 was not able to describe well the areas of the intersection. The main reason for this is that pxylene prefers the straight orientation (long axis of p-xylene parallel to the y direction, i.e., to the axis of the straight channels, marked as orientation C in Figure 10) when it finds itself in the intersection. As a result, it very rarely visits areas of the sinusoidal channel that are far from the red region within the intersection in Figure 10 and our sampling of these areas is not reliable. Repulsive walls constructed vertical to the sinusoidal channel path are not able to keep p-xylene in specific areas of the intersection, as we can see in Figure S5 of the Supporting Information. ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 42

When placed in these areas, p-xylene prefers to leave and visit areas of the sinusoidal channel which are not constrained by the repulsive walls. In addition, from the 3-D free energy profile we can see that constraining walls must be constructed with a large distance between them, so that they do not influence sampling in the intersection. In case they influence the sampling, the 1-D free energy fragments will not stitch perfectly together. After making these observations, we have tried to use other paths, which are illustrated with bold purple and bold red color in Figure S6 of the Supporting Information. The 1-D free energy fragments we obtained upon moving the repulsive walls along these paths were not able to be stitched together because they did not fit in some areas (see Figure S6 in the Supporting Information). If we examine Figure S7 of the supporting information, we can attribute this stitching problem to the existence of different orientational states. This means that we were not able to sample configuration space with precision, because p-xylene could be trapped in a free energy basin with low probability of escape.

Figure 10. 2-dimensional cross-section (at plane y=14.94Å) through the 3-dimensional free energy profile of p-xylene inside the sinusoidal channel for the DHFF force field and temperature T=300K. Red color shows the areas where the p-xylene center of mass spends most of its time. The colorbar is in RT units. In view of these problems, we chose to use a new path shown in red in Figure 10, instead of the old path shown in blue (actually including a spline to deal with the discontinuity in the intersection – see ref 27) in the same figure, and to deal with the issue of orientation within the intersection as discussed below. The difference between the old and the new paths is localized ACS Paragon Plus Environment

Page 17 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


just in the intersection area. This means that the 1-D p-xylene free energy profile can be compared directly with that of benzene for values of the ξ reaction coordinate along the sinusoidal axis from 3.2 to 10 Å. The new red path is shorter than the old blue one. This new red path is the reaction coordinate of the 1-D free energy profiles of p-xylene shown in Figure 11. We will elaborate on this Figure later.

Figure 11. i) Free energy profiles of p-xylene, in RT units, as a function of the reaction coordinate for 300K (blue), 405Κ (red), and 485Κ (green) along the sinusoidal channel of silicalite-1 for a crystal described by DHFF force field (a) and the SGVFF force field (b) and ii) free energy fragments of the area between the walls (where there is no influence of the bias potential of the wall), from which the 1-dimensional local free energy profile has been derived by the stitching process



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 42

As we mentioned, the probability of p-xylene to visit areas in the intersection with orientations different from the orientation it adopts inside the straight channels is small. Consequently, a new technique must be used to constrain p-xylene to these orientations, so that we can sample them sufficiently. To do this, we have invented orientationally constraining walls. First, these walls must keep p-xylene in specific orientations without influencing the dynamics, as was done in the case of the center of mass constraining walls we have already discussed. Secondly, these orientational walls must exert zero total force on p-xylene. If we use a bias potential V(y17), where y17=y7-y1, then a force will be exerted on the 1st and 7th atom only (see Figure 12) and only in the y direction (which seems to describe well the dividing surface of the orientational free energy-for the importance of the dividing surface, see Figure S1 of the supporting information). This force is described in detail by Eqs S2 and S3 of the supporting information.

Figure 12. Description of the orientation wall and the way we can obtain the 1-D orientation free energy profiles. Green color is used for the methyl united atoms. The repulsive potential V(y17) has the same form as the repulsive walls on the center of mass described in the work of Kolokathis et al.27 The y direction we chose for the orientation walls seems to describe well the orientational free energy dividing surface (see Figure S7 of the supporting information). Moving these orientational walls along the y17 coordinate, we obtain 1D orientational free energy profiles. To sample y17, we bin it in steps of 0.01 Å. In the projection of the center of mass of p-xylene and benzene, we have used steps of 0.05Å instead of 0.01Å. Because y17 can have values from -1.51Å to 1.51Å (very short compared with center of mass paths), we do not face noise issues using 0.01Å. As a result, this choice helps us improve our stitching accuracy. The 1-D orientational free energy profiles we get applying this procedure are shown in Figure 13. ACS Paragon Plus Environment

Page 19 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 13 betrays the existence of two orientational states in the intersection. The first state is the domain ξ[1.51 Å,0.90 Å)(0.90 Å, 1.51 Å] and we symbolize it as IS, while the second state is the domain ξ[0.90Å, 0.90 Å] and is symbolized as IZ. In the first state p-xylene adopts the orientation it likes when it is in the straight channel, depicted in Figure 7c, in Figure 8 and in Figure 10 (orientation C). This orientational state is much more favorable than the second state. In the second state, p-xylene prefers to have its ring in the x-z plane. Because this state is unfavorable and p-xylene rarely visits it, we decided to use orientation walls located in the saddle points of Figure 13. Simultaneously, we conducted molecular dynamics simulations with the use of center of mass repulsive walls for ξ belonging to [10.00Å, 13.00 Å]. Using this method we managed to obtain the 1-D free energy profile of Figure 11.

Figure 13. i) Profiles of the free energy of p-xylene, in RT units, as a function of the reaction coordinate at 300K (blue), 405Κ (red), and 485Κ (green) for a crystal described by the DHFF force field (a) and the SGVFF force field (b) and ii) free energy fragments of the area between the walls (where there is no influence of the bias potential of the wall), from which the 1dimensional local free energy profile has been derived by the stitching process. The reaction coordinate is the projection of a C-CH3 bond onto the y axis. ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 42

e) Rate constants, diffusion coefficients, and isosteric heat calculations From the 1-D free energy profiles (Figures 5, 6, 7, 9, 11 and 13 ) we can calculate the transition rate constants using eq 1 as soon as we define the states (areas between two maxima of the free energy). In eq 1, ( 1( i ) ,  2( i ) ) is the interval, along the reaction coordinate, corresponding to state i and  † is the reaction coordinate value at the free energy maximum separating states i and j. m is the mass of the diffusing molecule. In the cases we studied, a large number of states have been found. To be able to compare exactly the transition rate constants calculated here using the SGVFF force field with the rate constants calculated using the DHFF force field, we propose a reduction of the total number of macro-states to 4. These states are the Intersection state in which the sorbate has the orientation of the straight channel (IS), the Intersection state in which the sorbate has the orientation of the zig-zag channel (IZ), the Straight Channel state (S), and the Zigzag or sinusoidal channel state (Z). In the case of benzene, the total number of macro-states is 3, as we have only one state in the intersection (I). If the suggested reduction in the number of states is correct, it should not influence the predicted values of the diffusion coefficients. Also, this reduction has the advantage of simplifying the description of the systems we study without losing any important information.

 A( † )  exp     k BT 

1/2

TST

k i j

 kT   B   2 m 

2( i )

 A( )   (i ) exp   kBT  d 1



1/2

 kT   B   2 m 

 2( i )



†

 d 

(1)

1( i )

In the sketch of Figure 14 two cases are illustrated, which we encountered in our 1-D free energy profiles. In the first case, two states with low energy barriers between them (< 1RT) and four transition rate constants are reduced to one state and two transition rate constants (see Figure 14a). In the second case, two states with high energy barriers (states 1 and 3) divided by a state with low energy barrier (state 2) are reduced to two states with high energy barriers (see Figure 14b). In addition, the new rates produced after the reduction must also satisfy the microscopic reversibility condition. A proof that this is the case for the states defined here is given in the Supporting Information (eq S1). To our knowledge, this is the first time this method is used for rate constant calculation after a reduction in the number of states. In the first case, the transition rates are calculated for the values ξ1 tο ξ3 of the reaction coordinate for the state 1 of Figure 14a. In the second case, a special manipulation is needed. Because of the high free energy of state 2 of Figure 14b, the sorbate resides there with low probability. In case the sorbate enters state 2, it will move on quickly to another state, because the free energy barrier that must be overcome in order to escape from state 2 is low. The probability for a transition from state 2 via the ξ2 saddle point to state 1 to take place will be different from the probability for a transition from state 2 via the ξ3 saddle point to state 3 to take place. The two probabilities are described by eq 2. Pξi is the conditional probability that a transition out of state 2 will take place via the route of ξi and ρ is proportional to the probability density of the sorbate center of mass. In other words, a fraction of


Page 21 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the trajectories which lead to transition from state 1 to state 2, will manage to pass successfully to state 3. This, of course, has similarities with the transmission coefficient definition.70,71

P2 



   2

, P3 

2

3



, where i  e

3

   2



A1 D ( i ) RT

(2)

3

Figure 14. Two different methods of reduction in the number of states are illustrated. (a) two states with low energy barriers between them (< 1RT) and four transition rate constants are reduced to one state and two transition rate constants, (b) two states with high energy barriers (states 1 and 3) divided by a state with low energy barrier (state 2) are reduced to two states with high energy barriers As a result, the rate constant for the transition from state 1 of Figure 14b to state 3 of the same Figure will be described by eq 3.

k1TST 3

  1/2  k BT  TST   k12 P3     2 m   

    2    3  2   2  3  1  d 

(3)

Combining the 1-D free energy profiles of the Figures 4, 5, 6, 8, 10, 12 and the method of Figure 13, we get the macro-states which are listed in Table S8 of the Supporting Information. For these macro-states, using eq 3, we calculate the transition rate constants. Solving a linear system of



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 42

TST TST equations obeying microscopic reversibility  k i j Pi  k ji Pj  , we calculate the equilibrium

probability Pi for each macro-state i. To compute the isosteric heat of sorption, we place repulsive walls at the entrances to each state (maxima in the 1-D free energy profiles) to trap the sorbate in this state. Then, we write out the sorbate-silicalite-1 potential energy (VS-Z) in every timestep of an MD trajectory confined to the state. From the recorded values we compute a mean energy of interactions for each state. The isosteric heat is then obtained from eq 4.66 states

Qst   VS  Z

NVT

i  RT , or Qst    Pi VSstate Z

i 1

  2 1 L  d D  lim  min Α          





   

NVT

 RT

2    

(4)

(5)

For the diffusion coefficient calculation, we have used the MESoRReD eq 5 for a variety of ν, where ν determines the sequence of 2 unit cells arranged along the considered direction with periodic boundary conditions at the ends. We also use Ld to denote the length of the unit cell in the considered direction. Very small eigenvalues are difficult to calculate and, consequently, we are limited to specific trials of ν (usually for ν in the interval (1, 15)). The plateau value with respect to  is the diffusion coefficient we are looking for.28 To apply eq 5, we first assign identity numbers to each macro-state in the unit cell. Then we collect all rate constants pertaining to a unit cell in a matrix K1. For example, the rate k TST will be the element K1(j,i) of the K1 matrix. The i j

diagonal elements of this K1 matrix are K1 (i, i )   K1 ( j , i ) . For a specific direction x, y or z, a j i

KR1 matrix is constructed from K1. KR1 contains the rate constants of transition from a unit cell to the cell neighboring it on the right along the considered axis. All its remaining elements are zero. Similarly, a KL1 matrix is constructed containing only the rate constants of transitions between a unit cell and its neighbor to the left. Knowing K1, KR1, and KL1, one constructs the matrix   π   π  A ν  K1  1  cos  ν 1    K R1  K L1   i sin  ν 1   K R1  K L1  , with i being the imaginary unit.  2  2  

 

The quantity min Α appearing in eq 5 is the smallest among the absolute values of nonzero eigenvalues of matrix A ν , We collect all the results for rate constants, probabilities, isosteric heats and diffusion coefficients in Tables 1-4 for p-xylene-silicalite-1 (DHFF), and p-xylene-silicalite-1 (SGVFF) and in Tables S11-S14 of the supporting information for the benzene-silicalite-1 system. We have also used eq 6, proposed in ref 72, to calculate the expected p-xylene diffusion coefficient in the z direction of silicalite-1 (Dzz,t) based on the Dxx, Dyy values we calculated from eq 5. Then we compared these values with the Dzz values derived from eq 5, or directly from ACS Paragon Plus Environment

Page 23 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molecular dynamics simulations via mean square displacement versus time diagrams (see Tables 4a, 4b and Tables S14a, S14b of the Supporting Information). From this comparison we concluded that eq. 6 is satisfied to an excellent approximation by our simulation results for both benzene and p-xylene. This shows that the assumption of short particle memory, on which eq. 6 is based, is valid.72

c2 a2 b2 = + Dzz ,t Dxx Dyy

Table 1a. Rate constants, in s-1, of transitions of p-xylene in silicalite-1 (DHFF model) Transition From IS S Z

IZ

Temperature 300K

To

405K

485K

IZ

2.34×108

1.74×109

5.61×109

S

4.47×107

5.60×108

1.65×109

IZ

6.59×107

1.25×109

4.08×109

IZ (route a)

7.82×107

7.86×108

2.25×109

IZ (route b)

9.20×107

7.11×108

1.88×109

Z(route a)

4.20×109

1.19×1010

1.56×1010

Z(route b)

4.94×109

1.08×1010

1.31×1010

IS

1.87×1011

2.90×1011

2.81×1011

*

Routes a and b are defined in Figure 10


(6)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 42

Table 1b. Rate constants, in s-1, of transitions of p-xylene in silicalite-1 (SGVFF model) Transition From IS S Z

IZ

Temperature 300K

To

405K

485K

IZ

4.74×107

1.46×109

4.79×109

S

7.25×106

1.50×108

4.35×108

IZ

9.86×107

1.31×109

3.71×109

IZ(route a)

9.77×107

8.05×108

2.45×109

IZ(route b)

6.19×107

5.43×108

1.40×109

Z(route a)

1.14×109

4.29×109

8.91×109

Z(route b)

7.25×108

2.89×109

5.09×109

IS

7.42×1010

1.83×1011

2.32×1011

*

Routes a and b are defined in Figure 10

Table 2a. Residence probabilities of p-xylene in silicalite-1 (DHFF model) Probabilities

Temperature 300K

405K

485K

IS

0.5724

0.6485

0.6395

IZ

7.16×10-4

0.0039

0.0128

S

0.3884

0.2885

0.2589

Z

0.0385

0.0592

0.0889

Table 2b. Residence Probabilities of p-xylene in silicalite-1 (SGVFF model) Probabilities

Temperature 300K

405K

485K

IS

0.9245

0.8586

0.8246

IZ

5.9×10-4

0.0069

0.0170

S

0.0680

0.0980

0.0965

Z

0.0069

0.0365

0.0619


Page 25 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3a. Mean Potential Energies and isosteric heat of sorption of p-xylene in silicalite-1(DHFF model) Potential Energies

Temperature

kJ/mol

300K

405K

485K

IS

-69.02

-65.48

-62.66

IZ

-53.21

-49.89

-47.65

S

-76.11

-72.33

-69.98

Z

-70.55

-67.97

-65.88

Qst

74.31

70.92

68.69

Table 3b. Mean Potential Energies and isosteric heat of sorption of p-xylene in silicalite-1 (SGVFF model) Potential Energies kJ/mol

300K

Temperature 405K

485K

IS

-68.95

-65.70

-63.06

IZ

-53.86

-49.00

-47.56

S

-70.38

-68.07

-65.92

Z

-66.96

-64.73

-63.10

Qst

71.51

69.15

67.10

*

Qst=70kJ/mol for p-xylene in rigid silicalite-1 at 250K66

**

experimental values Qst  [64.4,80.0] kJ/mol66



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 42

Table 4a. Self-Diffusion coefficients of p-xylene in silicalite-1 (DHFF model) 2 -1

Self-Diffusivity (m s )

Temperature (MD-see Figure S9 of Supp.Info)

Temperature (TST) 300K

405K

485K

485K

Dxx

1.641×10-12

2.225×10-11

9.179×10-11

1.047×10-10

Dyy

1.271×10-11

1.800×10-10

5.242×10-10

4.398×10-10

Dzz

6.239×10-13

8.285×10-12

3.221×10-11

2.945×10-11

Dzz,t

6.509×10-13

8.868×10-12

3.500×10-11

3.792×10-11

D

4.991×10-12

7.020×10-11

2.160×10-10

1.913×10-10

Table 4b. Self-Diffusion coefficients of p-xylene in silicalite-1 (SGVFF model) Self-Diffusivity (m2s-1) 300K

Temperature 405K

485K

Dxx

2.643×10-13

1.194×10-11

5.558×10-11

Dyy

3.327×10-12

6.395×10-11

1.780×10-10

Dzz

1.072×10-13

4.369×10-12

1.823×10-11

Dzz,t

1.096×10-13

4.509×10-12

1.900×10-11

D

1.232×10-12

2.680×10-11

8.390×10-11

f) Investigation of entropic effects From the tables of the mean potential energies for each state (Tables 3 and S13) we can see that p-xylene has lower potential energy than benzene. We also have in mind that A  U  TS  K  V NVT  TS where U is the internal energy, K is the kinetic energy, V is the potential energy, T is the temperature and S is the entropy. In addition, we already know the Helmholtz free energy barrier ΔΑ from Figures 5, 7, 9. From now on, we define ΔΑ= Amax,IS - Amin,IS . Similarly, we define the barriers for the internal, the kinetic, the potential energies and the entropy. It is important to mention that ΔΚ=0 because the simulations are in the canonical ensemble and, consequently, (ΔU)NVT =(ΔV)NVT. From the previous analysis we see that we just need to calculate the average potential energy around the maximum (we remind the reader that we have already calculated the potential energy around the minimum in the intersection in Tables 3 and S13) if we want to calculate the entropic barriers for ACS Paragon Plus Environment

Page 27 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the transition from the intersection to the straight channel. To do this we use the repulsive walls to keep the center of mass of benzene and p-xylene in the interval ξ[12.60 Å,12.90 Å) of the straight channel. As a result, benzene and p-xylene move around the maximum of the 1-D free energy profile. At the same time, we write down every 1ps the potential energy between the sorbate and the crystal. The total time of the simulations is 5ns. From this sampling, we calculate the mean potential energy around the maximum of the 1-D free energy profile and finally we calculate the internal energy barrier. We performed the calculations for 300K only, because this is the common temperature between benzene and p-xylene in the simulations we mentioned before. We see these results in Figure 15 and Tables 5, 6. From Figure 15 it is clear that the entropy barriers for benzene are much higher than those for p-xylene, while the internal energy barrier contributes the most to the free energy barrier. We would expect Smax,Is

Diffusion of Aromatics in Silicalite-1: Experimental and Theoretical

Recommend Documents