Pressure- and Temperature-induced Monoclinic to Orthorhombic

experimental data. Thermodynamic considerations dictate that this is a second order phase transition in the Ehrenfest classification. Additionally, re...
3 downloads 8 Views 2MB Size
Subscriber access provided by - Access paid by the | UCSB Libraries

Article

Pressure- and Temperature-induced Monoclinic to Orthorhombic Phase Transition in Silicalite-1 Nikolaos Lempesis, Natalia Smatsi, Vlasis G. Mavrantzas, and Sotiris E Pratsinis J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b00400 • Publication Date (Web): 28 Feb 2018 Downloaded from http://pubs.acs.org on March 2, 2018

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 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 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 38 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

Pressure- and Temperature-induced Monoclinic to Orthorhombic Phase Transition in Silicalite-1

Nikolaos Lempesis,(1) Natalia Smatsi,(1) Vlasis G. Mavrantzas,(1),(2),* Sotiris E. Pratsinis(1) (1)

Particle Technology Laboratory, Institute of Process Engineering, Department of Mechanical

and Process Engineering, ETH Zürich, Sonneggstrasse 3, Zürich CH-8092, Switzerland (2)

Department of Chemical Engineering, University of Patras & FORTH-ICE/HT, Patras, Greece

*

Corresponding author: Vlasis G. Mavrantzas: [email protected]

Abstract: The thermal, mechanical, and volumetric behavior of silicalite-1, an all-silica MFI zeolite, is elucidated by atomistic simulations. A flexible force field was selected and validated from a set of force fields to capture the intramolecular interactions of the crystal lattice. This force field accounts for realistic bond, angle, and torsional interactions among atoms of the framework alongside with conventional Lennard-Jones and Coulomb interactions. By monitoring the behavior of silicalite-1 as a function of pressure and temperature, a fully reversible monoclinic to orthorhombic phase transition (polymorphism) was revealed in accordance with experimental data. Thermodynamic considerations dictate that this is a second order phase transition in the Ehrenfest classification. Additionally, reversible pressure-induced amorphization was captured by our model and was associated with the formation of linear zones of increased distortion running parallel to the straight and sinusoidal channels of this zeolite. Remarkably high isothermal compressibility (small bulk modulus) was calculated for orthorhombic silicalite1, in excellent agreement with experimental data, rendering silicalite-1 the most compressible zeolite known to date. The rigid unit mode (RUM) model was identified as the dominant structural mechanism for negative thermal expansion (NTE), typically observed over a wide temperature range in MFI zeolites. Better understanding of the monoclinic to orthorhombic phase transition, and molecular mechanisms associated with energy dissipation and NTE in zeolites provides control over the framework microstructure, allowing for enhanced molecular

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

sieving, tunable selectivity in separation processes, mechanical stability, and substantially amplified catalytic efficiency in petrochemical applications.

1. Introduction Zeolites are microporous crystalline materials composed of vertex-sharing TO4 tetrahedra, where T-sites are typically occupied by silicon or aluminium atoms.1 Cations reside within the zeolite framework ensuring electroneutrality by compensating for the excess negative charge caused by the substitution of silicon (Si+4) atoms by aluminium (Al+3) atoms. Other elements, such as zirconium,2 titanium,3 zinc,4 and germanium5 can also be found in the framework, attributing different properties and functionalities to the zeolite. Zeolites have been attracting the interest of the scientific and industrial community since their discovery in 1756.6 This is evident by the large number of industrially-developed crystal framework types (192), of the 232 frameworks that are known to date.1, 7 Owing to their unique porous structure and Brønsted acidity,8 zeolites constitute a remarkably versatile family of materials exhibiting a broad variety of engineering applications, including separation processes,9-10 molecular sieving and neutralization of heavy metals,11 cleansing of radioactive waste,12-14 dehydration of liquefied natural gas (LNG),15-17 water treatment,18-20 as well as catalytic cracking,21-22 and reforming.23-25 Through suitable selection of operating conditions, the microporous structure and associated sieving properties of zeolites may be finetuned and tailor-made. For example, a switch from positive, at small temperatures, to negative, at larger temperatures, thermal expansion in zeolites has been related to trapped water in the framework.26-28 By increasing the temperature, the trapped water is removed, allowing for contraction of the dehydrated zeolite upon further heating due to its intrinsic negative thermal expansion (NTE). This results in volume reduction and pore constriction, thereby increasing the separation selectivity of H2/CO2 mixture for CO2 capture and storage and H2 purification.10 Therefore, understanding the thermal behavior of zeolites is crucial for refining performance and effectiveness of applications that involve temperature variations. Zeolite Socony Mobil-5 (ZSM-5) zeolites belong to the structure-type code Mobil-Five (MFI), and are characterized by an Al-doped siliceous framework with Si/Al ratio ranging from 2 ACS Paragon Plus Environment

Page 3 of 38 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

10 to infinity; when this ratio exceeds 1000, the zeolite is typically known as silicalite-1.29 Silicalite-1 contains a pore system consisting of two types of interconnected channels. The first type of channels is a straight 10-membered ring channel running through the crystal parallel to the [010] direction (b-crystallographic axis), whereas the second type of channels is a sinusoidal 10-membered ring channel running parallel to the [100] direction (a-crystallographic axis). A tortuous pore path is observed in the [001] direction (c-crystallographic axis). ZSM-5 zeolites exhibit polymorphism, as they crystallize in a) an orthorhombic crystal lattice30 with space group Pnma, b) a monoclinic lattice31 with space group P21/n11, and c) an orthorhombic crystal lattice32 with space group P212121. In what follows, the first two will be referred to as ORTHO and MONO, respectively. The latter polymorph is observed predominantly in adsorption processes, when relatively large molecules (kinetic diameter ~ maximum pore size) were adsorbed in the framework,32 which resulted to a symmetry change from monoclinic (P21/n11) to orthorhombic (P212121). Since this work focuses on the structural changes of silicalite-1 induced by temperature and pressure, this polymorph won’t be considered here. A monoclinic (P21/n11) to orthorhombic (Pnma) phase transition (MOPT) has been reported in ZSM-5 zeolites.33-35 The monoclinic structure is thermodynamically stable below the transition point and the orthorhombic above.36 In general, MOPT in ZSM-5 zeolites can be instigated by a variety of external stimuli, such as a) temperature,37-38 b) pressure (or stress),39 c) concentration of defects40 (e.g. point defects created by substituting Si lattice atoms by Al),41 and d) the nature and concentration of sorbate loading.32 The temperature-induced MOPT is reversible,38 and occurs between 330 and 370 K at atmospheric pressure2,

37, 39, 41

with no

apparent hysteresis. The transition temperature depends, however, on the framework heteroatoms composition. The introduction of aluminium lowers the transition temperature such, that for relatively small Si/Al ratios (less than 100), the transition temperature falls below room temperature.42 Similarly, introduction of other trivalent dopants (such as, e.g. boron and/or iron) also decreases the transition temperature.43 In contrast, the transition temperature increases with increasing germanium content.43 On the other hand, the pressure-induced MOPT was reported to take place between 1 and 1.5 GPa at constant room temperature,39,

44

while there is still

ambiguity in the literature whether it is reversible44-45 or not39 upon decompression. 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

According to the Ehrenfest classification of phase transitions, polymorphism constitutes a second order phase transition.46-47 In that respect, the Gibbs free energy G and all its first order partial derivatives with respect to pressure and temperature, which relate to volume V and entropy S respectively, are continuous at the transition point. Additionally, according to this definition, at least one of the second order partial derivatives of the Gibbs free energy with respect to temperature and pressure, which relate to the isothermal compressibility βT, volumetric thermal expansion coefficient αV, and heat capacity at constant pressure Cp, exhibit a discontinuous “jump” near the transition point. Therefore, by monitoring the change of volume with pressure at constant temperature (expressed by βT), the change of volume with temperature at constant pressure (captured by αV), and the change of enthalpy with temperature at constant pressure (described by Cp), one can identify potential discontinuities that mark the transition point. MOPT is accompanied by structural changes within the crystal lattice, which, in turn, influence the overall behavior and properties of the material. For example, the monoclinic polymorph has the traits of a ferroelastic material, while the orthorhombic polymorph is paraelastic in nature.48 Of particular industrial relevance, MONO is more suitable for adsorption and diffusion processes49 which involve relatively large molecules, due to its larger maximum pore dimensions (5.85 Å) and more pronounced elliptical pore shape.31 ORTHO, on the other hand, owing to its smaller (5.3 Å) and more circular pore dimensions is more effective for molecular sieving, separation processes, and shape selective catalysis.31 In parallel to the pressure-induced MOPT, pressure-induced amorphization (PIA)50 in silicalite-1 was investigated here. Greaves et al. have shown in their seminal works51-52 that, upon compression, zeolites may transform successively into two different amorphous phases: a) a low-density amorphous (LDA) phase at intermediate pressures, and b) a high-density amorphous (HDA) phase at large pressures.53-54 The occurrence of more than one amorphous phases with the same chemical constitution but rather different densities and degree of order is called polyamorphism and has been observed in a variety of solid systems under compression.55 Fu and co-workers56 have shown, that silicalite-1 experiences PIA and assumes progressively two amorphous phases (LDA and HDA phases) upon compression. The formation of the LDA phase starts at around 1.5 GPa, while at around 3 GPa the LDA phase transforms gradually into 4 ACS Paragon Plus Environment

Page 5 of 38 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

the HDA phase. Up to this point, the process is considered reversible and the semicrystalline system can revert back to the fully crystalline state. Between 1.5 and 11 GPa, prevalent crystalline domains coexist with the LDA phase, that transforms progressively to HDA phase. Above 11 GPa, only the HDA phase exists and the PIA is deemed irreversible. This two-step amorphization process has been associated with the flexibility of zeolites described by the socalled “flexibility window”.57 Another industrially-significant characteristic of ZSM-5 zeolites is their complex, anomalous and highly anisotropic thermal behavior.58 Consequently, ZSM-5 zeolites exhibit positive thermal expansion (PTE) at low temperatures (300 – 425 K),59 whereas they contract upon further heating, exhibiting, thus, NTE.60 The transition from PTE to NTE with increasing temperature has been confirmed by numerous experimental works,61 however, the underlying mechanism of the phenomenon remains, in large part, still elusive. A broad variety of different explanations have been proposed to interpret it, including a) dealumination60 (i.e. expulsion of aluminium from the framework), b) dehydration,27 c) intrinsic framework mechanism,2 and d) MOPT.58 Most recent works,62-63 argue for the superiority of the rigid unit mode (RUM) model;64 an intrinsic framework mechanism activated by increasing temperature. According to the RUM model,65 [SiO4]4- tetrahedra are deemed as rigid units with very stiff Si-O bonds.66 Upon heating, tetrahedra do not interpenetrate one another, but rather rotate in a concerted way to accommodate unit cell volume reduction,67 as shown schematically in Figure 1. By focusing on heating up orthorhombic silicalite-1, one can preclude the influence of: a) MOPT, b) dealumination (no Al in our ORTHO model crystal), and c) dehydration (no water present in our ORTHO model), and therefore isolate and quantify solely the effect of the RUM model on the NTE of silicalite-1.

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

Figure 1: Schematic representation of the RUM model. Cyan rhombuses depict 2D projections of [SiO4]4- tetrahedra. Silicon atoms are at the center and oxygen atoms at the vertices of each tetrahedron. Unit cells at T0 and T1 > T0 are represented by the yellow and purple rectangular areas, respectively. Upon heating, the tetrahedra rotate by an angle θ = θ(T). Vectors R1 and R2 are also shown in red and blue, respectively (see text for details).

This report is organized as follows. Section 2 describes the systems and methods employed to simulate silicalite-1 and is divided into four subsections. Section 2.1 presents the details of the model construction and the force field selection process, whereas section 2.2 provides a detailed description of the structural analysis implemented for structure and force field validation purposes, section 2.3 discusses the methods invoked to capture temperature- and pressureinduced MOPT, as well as PIA by MD simulations, and section 2.4 describes the simulated heating of orthorhombic silicalite-1 for studying its thermal behavior. Section 3 reports and discusses the results of this simulation work. Force field and structural validation is provided in Section 3.1. The findings of the pressure- and temperature-induced MOPT are presented in sections 3.2 and 3.3, respectively. Section 3.4 illustrates the results related to the thermal behavior of silicalite-1 and the identification of the RUM model as the dominant intrinsic framework mechanism for NTE in pure siliceous zeolites. Finally, section 4 summarizes the basic conclusions of this work and provides proposals for future research.

2. Systems and Methods 2.1 Model construction and force field 6 ACS Paragon Plus Environment

Page 7 of 38 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

Here, two silicalite-1 model crystals were considered. Both models represented perfect, defect-free crystals. The first one was monoclinic (MONO), whereas the second was of orthorhombic symmetry (ORTHO). The unit cells of MONO and ORTHO were built by Scienomics MAPS,68 a software tool for building and analyzing realistic multiscale materials models. In that process, the fractional coordinates and space group symmetry of MONO31 and ORTHO30 were used as input. To create the initial atomistic configurations of MONO and ORTHO, the corresponding unit cells were replicated in the a, b and c crystallographic directions, for a total of (4 x 4 x 4) unit cells. This resulted in two simulation boxes of total size 8.043 x 7.952 x 5.347 nm for MONO and 8.028 x 7.968 x 5.368 nm for ORTHO, each one containing 18432 sites (12288 oxygen and 6144 silicon atoms). The lattice parameters, lattice symmetry, and space group of both forms of silicalite-1 considered here are summarized in Table 1. Two-dimensional projections of the two unit cells and simulation boxes along the [100] direction (parallel to the sinusoidal channels) are depicted in Figure 2. The subtle difference in the α angle between the two unit cells can be clearly seen in Figure 2.

Table 1: Unit cell details of the two forms of silicalite-1 considered in this work Parameter

ORTHO30

MONO31

a, nm

20.070

20.107

b, nm

19.920

19.879

c, nm

13.420

13.369

α, o

90.00

90.67

β, o

90.00

90.00

γ, o

90.00

90.00

Lattice system

Orthorhombic

Monoclinic

Space group

Pnma

P21/n11

7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

Figure 2: Two-dimensional projection of the unit cell of ORTHO (upper left) and MONO (upper right) viewed along the [100] direction. Simulation boxes for ORTHO (lower left) and MONO (lower right) created by (4 x 4 x 4) unit cell replication in 3D-space. Periodic boundary conditions were applied in all three directions of the simulation boxes. The outline of a unit cell is shown (blue box). Shaded cyan areas mark the cross section of indicative sinusoidal 10membered ring channels. Silicon atoms are depicted with black and oxygen atoms with red.

To select an accurate force field for the silicalite-1 framework, a very thorough investigation was conducted. In the vast majority of past simulation works on silicalite-1, mostly non-bonded framework-framework interactions were included in the employed force fields,69-76 neglecting thereby, realistic motions arising from stretching of bonds, bending of angles and torsions. In the few exceptional cases, where flexible force fields were used, either silicon atoms were not taken into consideration in calculating non-bonded interactions between framework atoms,77 or some bonded interactions (usually torsions) were neglected.78-79 All of these approaches, albeit simple in their implementation, did not reproduce accurately the volumetric behavior of the framework.80 Contrary to that, we tested a set of flexible force fields, taking into consideration all possible bonded and non-bonded contributions to the total system energy. Besides the 8 ACS Paragon Plus Environment

Page 9 of 38 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

conventional non-bonded Lennard-Jones and Coulomb interactions, bonded interactions related to bond-stretching, angle-bending, and torsions were considered. The force field selection for silicalite-1 proceeded based on two different criteria: a) the ability of the force field to capture accurately the volumetric behavior and structure of the crystal, and b) the existence of interaction parameters for Al heteroatoms that may be present in the framework. To the best of our knowledge, none of the past simulation works on silicalite-1 was able to reproduce the experimental framework density30 ρexp = 1.785 g/cm3. Instead, past works focused primarily on transport properties of sorbate molecules in the zeolite framework81 and on sorption processes.76, 82 As described in Section 1, most applications involving silicalite-1 are based on its unique structure comprising pores and channels of varying size and geometry. Since this work aspires to identify the molecular mechanisms that determine the structure of silicalite-1, and since the structure imparts different properties and functionalities to the zeolite itself, an accurate description of its volumetric behavior is of paramount importance. For the validation of a suitable force field, the ORTHO crystal, for which the experimental density was available,30 was transferred to LAMMPS,83 with the help of which molecular dynamics (MD) simulations were realized in the NPT ensemble using the Parrinello-Rahman method which allows for independent changes both in the dimensions and angles of the simulation cell in the course of the MD simulation.84-85 The system volume changes freely and at the end of the simulation it assumes the equilibrium value dictated by the used force field. Totally seven widely used and thoroughly validated flexible force fields (Amber,86 Dreiding,87 Pcff,88 Cvff,89 Uff,90 Hill and Sauer,91 improved Hill and Sauer92) were used to fully equilibrate ORTHO with a time step of 1 fs. A cut-off distance of 10 Å was used for both Lennard-Jones and Coulomb interactions, as well as a tail correction. For each tested force field, one MD simulation ran for a total of 5 ns, the last 4 ns of which were sampled every 500 fs to calculate ensemble averages. Owing to the relatively quick equilibration in crystals, due to the dominant relaxation modes related to fast vibrational motions over thermal fluctuations, longer runs were not necessary. During the course of the MD simulations, the density of the system was monitored. The deterministic Nosé-Hoover thermostat and barostat were used to maintain isothermal and isobaric conditions,93-94 with time constants of 10 and 350 fs, respectively. The Velocity Verlet method95 was used to integrate the equations of motion. After completion of the 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

MD simulations, the average equilibrium density was calculated for every force field. As it will become apparent in section 3.1, the improved Hill and Sauer force field92 was found to be superior in capturing the volumetric behavior and overall structure of silicalite-1. For reasons of completeness, and given the importance of an accurate force field, tables containing the parameters and functional forms of all examined force fields are summarized in the Supporting Information. A point to note in these tables is that parameter B of the 9-6 Lennard-Jones interactions in the improved Hill and Sauer force field92 is zero for both types of atoms present in the framework, thus nullifying the attractive component of the Lennard-Jones interactions. However, this doesn’t make the improved Hill and Sauer force field solely repulsive, as attractive non-bonded contributions to the total energy arise from electrostatic (Coulomb) interactions between dissimilar atoms.

2.2 Structural analysis Having thus constructed an atomistically detailed model for ORTHO and MONO, and having selected an accurate force field with flexible contributions to the total system energy, the X-ray diffraction pattern of silicalite-1 was calculated. Due to software limitations, only the Xray diffraction pattern of ORTHO was extracted. This calculation served as an extra level of validation of the structure of the built model crystal, but also of the selected force field. For that purpose, the static structure factor S(q) of ORTHO was calculated with LAMMPS and compared to available experimental X-ray diffraction patterns for orthorhombic silicalite-1.7 The S(q) is a measure of how a material scatters incident radiation. The relation between S(q) and X-ray diffraction pattern has been described extensively elsewhere.96 For the calculation of S(q), the equilibrated structure of ORTHO was used as the initial configuration in a 1 ns long MD simulation in the NVT ensemble. For the calculation of S(q), the canonical ensemble (NVT) was chosen over the isobaric (NPT) to maintain the volume constant at its equilibrium value. Volume fluctuations influence the structure and, thus, the computed X-ray diffraction pattern. All other specifications for the MD run were identical to those mentioned in section 2.1. The Cartesian coordinates of all atoms of the crystal were stored every 0.1 ns along the aforementioned MD trajectory. An X-ray pattern was computed for every sampled configuration. The final X-ray diffraction pattern was obtained by averaging over all calculated diffraction patterns along the 10 ACS Paragon Plus Environment

Page 11 of 38 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

MD trajectory. The wavelength of the simulated incident beam used in this analysis was 1.54056 Å.

2.3 Monoclinic-to-orthorhombic phase transition in silicalite-1 Molecular dynamics simulations were implemented to study the MOPT of silicalite-1. The MOPT is accompanied by structural rearrangement within the crystal lattice of silicalite-1, changing thereby the average size, cross area, and shape of the pores and channels of the framework. Since this transition has immense consequences on most industrially-significant applications of zeolites, detailed MD simulations were employed to investigate the conditions that trigger it. Among the many external stimuli that can instigate the MOPT (see section 1), the effect of temperature and pressure on the structure of silicalite-1 were deemed as more meaningful and relevant to most practical applications. It is worth noting that the majority of engineering applications involving zeolites take place over a wide temperature and pressure range. Two series of MD simulations were realized. In each case, the equilibrated at ambient conditions MONO configuration was the starting point. In the first set of simulations, the effect of pressure on the structure of silicalite-1 was investigated by applying compressive stress in all three principal directions of the simulation box. We considered stresses from 0.01 to 5 GPa in increments of 0.25 GPa. This loading range is realistic in terms of zeolite synthesis39 and applications.97 The stress interval between 1.25 and 1.75 GPa, where the MOPT is expected to take effect at ambient temperature44, was sampled more densely in increments of 0.125 GPa. To examine the reversibility of the effect of pressure on the structure of silicalite-1, the aforementioned process was reversed in relevant decompression simulations starting from the 5 GPa configuration. For each stress, a 3 ns long MD simulation in the NPT ensemble was used at constant room temperature. The last 2 ns of every run were sampled every 0.5 ps for building ensemble averages. All simulations were run in triplicate with all other simulation details identical to those reported in section 2.1. In the second set of simulations, the effect of temperature on the structure of silicalite-1 was simulated by heating up MONO from 280 K to 400 K in increments of 10 K at constant 11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

atmospheric pressure. To check the reversibility of the effect of temperature on the structure of silicalite-1, the heating simulations were inverted by using the equilibrated configuration at 400 K as the starting point of a series of cooling simulations. Here also, triplicate MD simulations were realized in the NPT ensemble at each temperature. Each simulation ran for 5 ns, the last 3 ns of which were sampled every 0.2 ps for property estimation. To quantify the effect of PIA in silicalite-1, the Lindemann index was used.98 The Lindemann index is a measure of the averaged bond length fluctuations brought about by external stimuli99 such as temperature or pressure variations, and is used as a simple proxy for pressure-induced amorphization. In the solid state, the only motion that atoms execute is practically small vibrations around their equilibrium lattice points. As the material is pushed away from this state of equilibrium (by, e.g., compression), each atom is provided with energy that allows for elevated vibrations around the prescribed equilibrium lattice points. When these fluctuations exceed a certain threshold value, the structure should not be considered perfectly ordered anymore. The Lindemann index is particularly useful for capturing temperature-induced amorphization during melting of systems containing nanoparticles100, however, in this work, we used it to quantify the PIA of silicalite-1. For a system of N atoms, the system-averaged Lindemann index q is given by equation 1, where qi is the local (atomic) Lindemann index for the ith atom and rij is the distance between the ith and jth atoms. Angular brackets denote thermal averaging at temperature T.

q

1 N

q

(1)

i

i

1 qi   N -1 ji

rij2

T

rij

 rij

2 T

(2)

T

2.4 Negative thermal expansion (NTE) in silicalite-1: validating the RUM model Next, an NPT MD simulation with simulated heating was used, to investigate the thermal behavior of silicalite-1 over a broad temperature range. In order to avoid the possible influence of the temperature-induced MOPT on the overall thermal behavior of silicalite-1 (see section 1), 12 ACS Paragon Plus Environment

Page 13 of 38 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

the heating simulation was initiated from the equilibrated ORTHO configuration. The ORTHO was then heated from T0 = 298.15 K to T1 = 800 K over a 15 ns long MD simulation at constant atmospheric pressure. Also here, all remaining simulation details were identical to those mentioned in section 2.1. The Cartesian coordinates of all atoms were stored every 0.25 ps. To verify the efficiency of the RUM model in explaining the NTE in silicalite-1, two vectors were defined. The first vector R1 connected two neighboring silicon atoms, each one situated at the center of adjacent, vertex-sharing tetrahedra of the framework. Vector R2 was defined as the vector connecting the center of every tetrahedron, occupied by a silicon atom, with every one of its vertices, occupied by oxygen atoms. The two vectors R1 and R2 are depicted in Figure 1 with red and blue arrows, respectively. The average distance between neighboring silicon atoms, expressed by the magnitude of R1, was tracked during the heating simulation. Rotation of tetrahedra upon heating is expected to induce a reduction of the average distance between neighboring silicon atoms, as shown in Figure 1. Additionally, the normalized time autocorrelation function (ACF) of vector R2 corresponds to the cosine of the rotation angle θ (Figure 1). The ACF of vector R2 is calculated by equation 3:

ACF 

R2  t   R2  0

R2  0  R2  0

 cos θ

(3)

In eq 3, vectors R2 were transformed to unit vectors, whereas angular brackets denote averaging over all possible vectors R2 in the system and over all sampled times t. Angle θ depends on temperature and is expected to increase with increasing temperature owing to more intense thermal motion of atoms at elevated temperatures. However, it cannot exceed 90o as that would correspond to a perfect stacking of tetrahedra on top of one another (Figure 1). Thus, rotation of tetrahedra upon heating by θ (0ο ≤ θ < 90ο) would translate into a partial decay of the ACF with time. In this simulation, time and temperature are equivalent, as temperature increased monotonically from T0 = 298.15 K to T1 = 800 K over the 15 ns long heating simulation.

3. Results and Discussion 3.1 Force field validation and structural analysis 13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

As mentioned already, as first step of validating a force field, the volumetric behavior of silicalite-1 was examined. To this end, MD simulations were employed in the NPT ensemble to equilibrate the structure of silicalite-1 by using various flexible force fields. For each simulation, the equilibrium density of silicalite-1 was monitored and compared to the experimental value30 of ρexp = 1.785 g/cm3. Table 2 shows a comparison of the calculated equilibrium densities of silicalite-1 for all considered force fields. The force field benchmarking showed that the very recent modification92 of the original, well-known Hill and Sauer force field91 captured the density of silicalite-1 with remarkable accuracy (standard error < 1.5%). The subtle changes introduced by Boulfelfel et al.92 predicated on reestablishing the equilibrium values for the Si-OSi and O-Si-O angles. Thus, in what follows, the improved Hill and Sauer force field derived by Boulfelfel et al.92 was used for all MD simulations. At this point a clarification is in order. The implementation of the improved Hill and Sauer force field in LAMMPS proceeded by modifying the LAMMPS source code such that the 9-6 Lennard-Jones potential used in the Hill and Sauer force field is represented correctly.

Table 2: Equilibrium density of silicalite-1 estimated by different force fields in NPT MD simulations at T = 300 K and P = 1 atm. The standard error of the simulated value, its deviation from the target value and the ability of each force field to account for aluminum heteroatoms are also shown 30

Exp.

Density [g/cm3] St. error [g/cm3] Deviation Aluminum parameters

Amber

86

87

Dreiding

Pcff

88

Cvff

89

Uff

90

Hill and Sauer91

Improved Hill and Sauer92

1.785

1.838

2.545

1.670

1.876

1.884

1.677

1.762

-

0.040

0.052

0.054

0.056

0.048

0.046

0.034

-

2.97%

42.58%

6.44%

5.10%

5.55%

6.05%

1.29%

-

NO

YES

YES

NO

NO

YES

YES

As additional validation of the selected force field, as well as of the structure of the built model crystal, the X-ray diffraction pattern of ORTHO was calculated and compared with 14 ACS Paragon Plus Environment

Page 15 of 38 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

available X-ray diffraction patterns of orthorhombic silicalite-1 powder taken from the International Zeolite Association.7 In accordance with similar comparisons for crystalline polymers101 and coalescing bimetallic nanoparticles102, Figure 3 shows that the agreement between the MD-calculated diffraction pattern and its experimental counterpart is very good. The characteristic peaks of the X-ray diffraction pattern,103 marked with dots on top of the graphs, are captured successfully. To the best of our knowledge, this is the first time that the X-ray diffraction pattern, and thus the structure of silicalite was calculated via simulations. The very good agreement between simulation result and experiment further validates the ability of the selected force field to capture accurately any structural changes in the framework of silicalite-1. Additionally, this agreement ratifies that our model crystal represents the microstructure of silicalite-1 realistically. Since the spatial extent of the framework is dictated by the intramolecular interactions among atoms, capturing the equilibrium density and microstructure of silicalite-1 translates into an accurate representation of its actual pore size and channel geometry by our model crystal.

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

Figure 3: Calculated (black line) and experimentally measured7 (red line) X-ray diffraction pattern of orthorhombic silicalite-1. For clarity, both patterns are normalized. The characteristic peaks of silicalite-1 are marked with blue dots.103

3.2 Pressure-induced MOPT and PIA in silicalite-1 The response of silicalite-1 to pressure was assessed by MD simulations. The starting structure was the MONO configuration which was subjected to a set of pressures ranging from 0.01 to 5 GPa. The process was then inverted in a series of decompression runs. During each of these runs, the crystal’s volume was monitored. Figure 4 shows the calculated change of unit cell volume with pressure at room temperature. This change is quantified by the isothermal compressibility βT, or equivalently, by the first order partial derivative of volume with respect to pressure, the latter corresponding to the slope of V=V(P) of Figure 4. A change in slope of V=V(P), would correspond to a discontinuous change of βT. Figure 4 shows, that there are two 16 ACS Paragon Plus Environment

Page 17 of 38 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

distinct regions in the dependence of volume V on pressure P, characterized by different slopes of the V=V(P) curve. These regions are marked with different shading colors in Figure 4. For reasons of clarity, error bars have been omitted from Figure 4; the standard error in all cases is on the order of 2% of the mean. The change in slope of the V=V(P) curve occurs at a transition pressure Ptr ≈ 1.4 GPa. This observation constitutes a first indication that a second order phase transition takes place at around Ptr ≈ 1.4 GPa. The bulk modulus K of the ORTHO model crystal, was computed from the slope of V=V(P) after the phase transition point as the inverse isothermal compressibility K = βT-1. Our result (K = 18.5 ± 0.4 GPa) is consistent with experimental data of Quartieri et al.39 (K = 18.2 ± 0.2 GPa) and Haines et al.44 (K = 18.8 ± 0.5 GPa) for silicalite-1 synthesized in fluoride media and compressed with “non-penetrating” media, despite the existence of defects in the latter specimens caused by the presence of hydroxyl groups.40 Given that typical values of zeolite bulk moduli range from about 18 to 72 GPa,97 silicalite-1, exhibiting a remarkably low bulk modulus (K = 18.5 GPa), can be classified as the most compressible zeolite structure.39 This result suggests insurmountable mechanical stability of silicalite-1, which in conjunction with ongoing research efforts aiming at improved thermal and chemical stability104 can lead to extremely stable thin zeolite membranes, with immense impact on separation performance,105 and catalytic activity.106

Figure 4: Change of unit cell volume of the monoclinic configuration (MONO) as a function of pressure during compression (solid line) and decompression (broken line) at T = 300 K. Different shading colors are used to designate different slopes of V=V(P). 17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

To examine in more detail the structural change of the MONO configuration with pressure, as well as the associated possibility of a pressure-induced MOPT, we have monitored the change of lattice angle α of the monoclinic unit cell as a function of pressure. As mentioned in section 2.2, angle α is about 90.67o in the monoclinic configuration. Figure 5 shows that the value of α drops with increasing pressure, reaching that of the orthorhombic lattice (α = 90ο) at about Ptr = 1.4 GPa and remaining practically unaltered thereafter, retaining thus the orthorhombic symmetry. It is worth noticing that the deviation from the mean value of α becomes larger with increasing pressure, demonstrating that, as pressure increases, some unit cells of the simulation box experience larger distortion away from their equilibrium structure. Nevertheless, these distortions in the crystal lattice, as large as they might be, seem to cancel out resulting in an average value of practically unaltered α. By combining the results of figures 4 and 5, one can conclude that a pressure-induced symmetry change from monoclinic to orthorhombic takes place at Ptr ≈ 1.4 GPa in silicalite-1. This simulation result compares favorably to a series of experimental observations reporting a pressure-induced symmetry change in silicalite-1 between 1 – 1.5 GPa.39, 44 Nonetheless, it is still debatable in the literature whether this transition is reversible. However, by looking at our decompression results illustrated in both Figures 4 and 5, we conclude that the pressure-induced MOPT is reversible, within the uncertainties of the simulations.

Figure 5: Change of the unit cell angle α of the monoclinic configuration (MONO) as a function of pressure during compression (solid line) and decompression (broken line) at T = 300 K. 18 ACS Paragon Plus Environment

Page 19 of 38 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

Figure 6 shows the system-averaged Lindemann index as a function of imposed pressure. Three different regions can be distinguished in Figure 6; a) a low pressure region (0 – 1.5 GPa), where the Lindemann index remains practically constant, b) a narrow transitional region (1.5 – 2 GPa), where the Lindemann index increases rapidly, marking the onset of some kind of structural change with associated partial loss of order within the framework, followed by c) a broader region where the Lindemann index levels off reaching some kind of a plateau. The first, low-pressure region (P < 1.5 GPa), describes the equilibrium crystalline structure of silicalite-1, where atoms vibrate about their equilibrium lattice points in the monoclinic framework. No significant change in the Lindemann index is observed in this low-pressure region. As pressure increases, the MOPT occurs at around Ptr ≈ 1.4 GPa (Figures 4-5). It is worth noting, that the MOPT has no significant impact on the Lindeman index as the system retains its crystallinity and all atoms continue to assume their prescribed lattice points. Therefore, no loss of order is brought about by the MOPT. Figure 6 suggests that shortly after the MOPT some moderate amorphization takes effect at intermediate pressures (1.5 < P < 2 GPa). In fact, it has been demonstrated in the past that a LDA phase coexisted with the crystalline phase at intermediate pressures.107 This is validated by our results in Figure 6, where the MOPT is followed by a mild amorphization process reflected by a sudden increase of the system-averaged Lindemann index q. The third region of Figure 6 (P > 2 GPa) is related to the establishment of a low-density amorphous (LDA) phase. Interestingly, the crystalline phase remains still dominant, as manifested by the average angle α that maintains the orthorhombic value (α = 90o) throughout this pressure range (Figure 5). In fact, the increasing error bars of angle α as pressure grows, may be attributed, at least partially, to the existence of some LDA zones within the silicalite-1 framework. At this stage, the system can recrystallize fully back to its perfect crystalline structure upon pressure release, as confirmed by our decompression simulations. Additionally, no apparent hysteresis was observed during the decompression process, manifesting that no dissipative mechanisms came into effect at the considered pressure range.

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

Figure 6: System-averaged Lindemann index of silicalite-1 as a function of pressure at T = 300 K during a) compression (solid line) and b) decompression (broken line). Different shading colors are used to mark the three distinct regions (see text for details).

By examining the behavior of the local (atomic) Lindemann index with pressure, we visualized the onset of amorphization between 1.5 and 2 GPa. For larger pressures (P > 2 GPa), the amorphization process was associated with the formation of linear zones of increased qi parallel and perpendicular to the straight and sinusoidal channels. This becomes apparent by looking at the heat maps provided in Figures 7 and 8. A number of interesting features can be observed there. At relatively low pressures (P ≤ 1.5 GPa), Figures 7a and 8a exhibit practically a homogeneous distribution of atomic Lindemann indices, indicating that all atoms of the crystal framework accommodate the imposed pressure in a similar way. At intermediate pressures (1.5 < P ≤ 2 GPa), Figures 7b and 8b show the onset of “hot” zones formation, that are parallel to the [010] (straight channels) and [100] direction (sinusoidal channels), respectively. At larger pressures (P > 2 GPa), these linear zones spread throughout the material and become even more pronounced (Figures 7c-d and 8c-d), while they are complemented by the formation of additional linear “hot zones” in the [001] direction, that is, perpendicular to both straight and sinusoidal channels (not shown). These linear “hot zones” are associated with increased distortion of atoms away from their equilibrium lattice points in response to imposed pressure and can be rationalized as areas of either increased stress concentration or increased stress dissipation. In either case, this conclusion constitutes a first indication of an atomic mechanism in response to 20 ACS Paragon Plus Environment

Page 21 of 38 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

increased compressive pressure and may be deemed as a precursor mechanism for the onset of plastic deformation and eventually material failure. A more thorough analysis on the mechanics of silicalite-1 involving calculations of stress and free volume spatial distributions is, however, needed.

Figure 7: 2D projections viewed in the [100] direction (along the sinusoidal channels) of the atomic Lindemann index calculated via eq 2 at a) 1.5, b) 1.75, c) 2.25 and d) 4.75 GPa.

21 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

Figure 8: 2D projections viewed in the [010] direction (along the straight channels) of the atomic Lindemann index calculated via eq 2 at a) 1.5, b) 1.75, c) 2.25 and d) 4.75 GPa.

3.3 Temperature-induced MOPT in silicalite-1 The change of unit cell volume and enthalpy of the MONO configuration with temperature are depicted in Figure 9. The change of volume with temperature at constant pressure is quantified by means of the volumetric isothermal compressibility αV. Similarly, the change of enthalpy with temperature at constant pressure is described by the heat capacity Cp. A discontinuous change of either of these quantities relates to a second order phase transition. By contrast to the clear discontinuous change of βT at Ptr (Figure 4), Figure 9 shows no clear evidence of such a change in the slope of either V=V(T) or H=H(T), which change rather monotonically with temperature. This observation, however, is not deterrent of a second order phase transition. As mentioned in section 1, a discontinuous change of at least one of the second order partial derivatives of G with respect to T and P suffices for the identification of a second order phase transition in the Ehrenfest classification.

22 ACS Paragon Plus Environment

Page 23 of 38 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

Figure 9: Change of unit cell volume V (blue) and enthalpy H (red) of the monoclinic configuration (MONO) as a function of temperature at P = 1 atm during heating (solid line) and cooling (broken line). The slopes of V=V(T) and H=H(T) relate to αV and Cp, respectively. For reasons of clarity all error bars have been omitted; the standard error in all cases was less than 3% of the mean.

Figure 9 reveals the negative slope of V=V(T) that relates to a NTE coefficient αV, typical in MFI zeolites. The molecular origins of this phenomenon are examined and commented upon later on in this work. By calculating the slope of V=V(T) in Figure 9, the volumetric thermal expansion coefficient αV of silicalite-1 was calculated αV = -8.81 X 10-6 K-1. That compares favorably to experimental data: (-7.64 X 10-6 K-1),2 (-15 X 10-6 K-1),58 (-14.5 X 10-6 K-1).62 It is worth mentioning that values for experimental volume thermal expansion coefficients fluctuate strongly. This strong deviation is attributed to the existence of crystallographic defects in the framework related to impurities, traces of heteroatoms, as well as trapped water molecules in the framework pores.27 The lack of a positive thermal expansion in Figure 9 (positive slope of V=V(T)) is justified by the absence of water molecules and hydroxyl groups from our model crystal. The temperature-induced MOPT in silicalite-1 was validated by monitoring the change of the obtuse lattice angle α with temperature. Figure 10 shows how angle α changes with temperature during the heating and cooling simulations. It changed with increasing temperature from the monoclinic value α ≈ 90.67o to the orthorhombic α ≈ 90o. The transition temperature Ttr is identified as the inflection point of a hyperbolic tangential fit108 (not shown) to the simulation curves of Figure 10. Accordingly, the phase transition takes place at around Ttr = 340 K, which is 23 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

in agreement with experimental observations of the transition temperature ranging from 320 to 370 K.2,

37, 39, 41

As Figures 9 and 10 show, our simulations validated the reversibility of the

temperature-induced MOPT. Capturing the temperature-induced MOPT by monitoring angle α, while, at the same time, no distinct slope change of V=V(T) and/or H=H(T) was observed, proclaims that the effect of temperature on the structure of silicalite-1 is very subtle and mild in nature and, in any case, less aggressive than the effect of pressure discussed in section 3.2. We believe that this is probably the reason why temperature-induced phase transition in silicalite-1 is undeniably reported to be reversible,38 as opposed to the pressure-induced phase transition, whose reversibility was up to date still debatable.39, 44

Figure 10: Change of the unit cell angle α of the monoclinic crystal (MONO) as a function of temperature at P = 1 atm during heating (solid line) and cooling (broken line).

3.4 NTE in silicalite-1: Validation of the RUM model The thermal behavior of silicalite-1 was monitored with the help of a MD simulation, during which orthorhombic silicalite-1 was heated up from T0 = 300 K to T1 = 800 K. On the course of this simulation, we monitored the change of interatomic distances described by vectors R1 and R2 (Figure 1 and discussion in section 2.4). Figure 11a shows the decay of the normalized time autocorrelation function of vector R2 calculated by virtue of eq 3. The time ACF of R2 decorrelates partially with increasing temperature, leveling off at around 0.4. Since the ACF corresponds to the average cosine of the tetrahedra rotation angle θ, the observed behavior illustrated in Figure 11a denotes a relative decrease of from its initial value 1 (θ 24 ACS Paragon Plus Environment

Page 25 of 38 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

= 0ο) to a reduced value of about 0.4 (θ ≈ 68ο). In other words, we found that upon heating up orthorhombic silicalite-1 from T0 = 300 K to T1 = 800 K the [SiO4]4- tetrahedra rotated on average by an angle of about 68o, thus causing a reduction of the unit cell volume, confirming the schematic rationalization of Figure 1.

(a)

(b) Figure 11: (a) Normalized time autocorrelation function (ACF) of vector R2, and (b) the magnitude of vector R1 describing the distance between neighboring Si atoms over the 15 nslong MD heating simulation from T0 = 300 K to T1 = 800 K.

To further validate the rotation of tetrahedra by virtue of the RUM model, we monitored the magnitude of vector R1, connecting two neighboring Si atoms. Figure 11b shows the average magnitude of vector R1 on the course of the heating MD simulation. It is obvious, that the average distance between adjacent Si atoms decreases monotonically, indicating that neighboring Si atoms are approaching one another as temperature increases. It is worth mentioning that the fluctuations of the magnitude of vector R1 become more intense as temperature increases, 25 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38

indicating stronger thermal vibrations at elevated temperatures. Figure 11, in its entirety, provides compelling evidence of the validity of the RUM model in describing the NTE in silicalite-1. Since our model crystal is free of water and hydroxyl groups and already of orthorhombic symmetry (no MOPT), through this analysis, we have been able to identify exclusively the rotation of rigid [SiO4]4- tetrahedra as the sole intrinsic framework mechanism explaining NTE over a wide temperature range (300 – 800 K). The rotation of [SiO4]4- tetrahedra along the course of the heating simulation is depicted in Figure 12 for a subdomain of the 4x4x4 simulation box. In Figure 12, for clarity and brevity, bonding of the depicted subdomain with all neighboring atoms has been omitted. It can be seen clearly, that upon heating, neighboring tetrahedra rotate in a concerted fashion, as originally assumed by the RUM model (cf. Figure 1).

a

b

c

d

e

f

Figure 12: Rotation of vertex-sharing [SiO4]4- tetrahedra on the course of the heating simulation. Snapshots of a simulation box subdomain are shown at a) t = 0, b) t = 1, c) t = 2, d) t = 6, e) t = 9 and f) t = 10 ns. Facets of two neighboring vertex-sharing tetrahedra are outlined with yellow line, whereas vector R2 is shown in all parts. The color code used is consistent with that of Figures 1 and 2.

26 ACS Paragon Plus Environment

Page 27 of 38 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

4.

Conclusions Given the great practical utility and industrial significance of zeolite frameworks, operative

over a very wide range of temperatures and pressures, the thermal and mechanical behavior of silicalite-1 was investigated by atomistic simulations. For the first time, a flexible force field with realistic three-body and four-body interactions has been shown to capture accurately the volumetric behavior and equilibrium structure of silicalite-1. The equilibrium density and calculated X-ray diffraction patterns compared favorably to experimental measurements. The MOPT was thermodynamically identified as a second order phase transition (polymorphism) according to the Ehrenfest classification and was examined under variations of pressure and temperature. The pressure-induced MOPT was captured accurately and found to be reversible, as confirmed upon reversal of the compressive pressure. A discontinuous change of the isothermal compressibility βT, as well as lattice angle change from the monoclinic (90.67o) to the orthorhombic value (90o) was observed at Ptr ≈ 1.4 GPa, in agreement with experimental data reporting a pressure-induced symmetry change of silicalite-1 between 1 – 1.5 GPa. A remarkably high isothermal compressibility βT = 18.5 GPa-1 was calculated for orthorhombic silicalite-1 and found to be in excellent agreement with data in the experimental record, according to which silicalite-1 is the most compressible zeolite known to date. This increased mechanical stability suggests that silicalite-1 is hardly distorted upon sorption. Reversible pressure-induced amorphization was confirmed to start taking effect immediately after the pressure-induced MOPT, becoming however more pronounced with increasing pressure. For pressures larger than 2 GPa, a LDA phase was observed to coexist with the prevalent crystalline phase. Visualization of 2D facets of the crystal under compression, revealed the formation of linear zones of highly distorted atoms away from their equilibrium positions. These linear zones were parallel (at intermediate loadings) and perpendicular (at higher loadings) to the straight and sinusoidal channels providing preliminary indications of intrinsic framework mechanisms that become operative to accommodate dissipation of excess energy. Temperature-induced MOPT was identified to take place at Ttr ≈ 340 K, in agreement with experimental data (320-370 K). Our simulations provided a negative volumetric thermal expansion coefficient αV, which compared favorably to experimental values. The lack of a 27 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 28 of 38

positive thermal expansion coefficient at relatively low temperatures was attributed to the absence of water molecules and hydroxyl groups, typically present in experimental studies. By starting our heating simulations from orthorhombic silicalite-1, we were able to preclude contributions to the NTE coming from the MOPT, dealumination, and dehydration of the zeolite, thereby identifying and isolating the RUM model as the sole intrinsic framework mechanism for NTE in zeolites. In summary, our model crystal described very accurately the volumetric, structural, thermal and compressive behavior of a pure-siliceous zeolite, setting the foundations for further research on simulating industrially relevant applications (such as adsorption, diffusion and separation processes) involving pure or doped zeolite frameworks.

5.

Supporting Information

Supporting Information: Functional forms and parameters of examined force fields

6.

Acknowledgement

This work was supported by computational time granted from the Greek Research & Technology Network (GRNET) in the National HPC facility ARIS under the project name COACERVATE (pr004011). NL and VGM acknowledge financial support through the European Project FORCE (=Formulations and Computational Engineering) funded by the European Commission (H2020NMBP-2016-two-stage), Project No 71027.

7.

Abbreviations

ACF: autocorrelation function, HDA: high-density amorphous, IZA: International Zeolite Association, LDA: low-density amorphous, MD: molecular dynamics, MONO: monoclinic ZSM-5, MOPT: monoclinic-to-orthorhombic phase transition, NTE: negative thermal expansion, ORTHO: orthorhombic ZSM-5, PIA: pressure-induced amorphization, PTE: positive thermal expansion, RUM: rigid unit mode.

28 ACS Paragon Plus Environment

Page 29 of 38 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

8.

References

(1) Auerbach, S. M.; Carrado, K. A.; Dutta, P. K., Handbook of Zeolite Science and Technology; Taylor & Francis, 2003. (2) Bhange, D. S.; Ramaswamy, V., Negative Thermal Expansion in Silicalite-1 and Zirconium Silicalite-1 Having Mfi Structure. Mater. Res. Bull. 2006, 41, 1392-1402. (3) van der Waal, J. C.; Rigutto, M. S.; van Bekkum, H., Zeolite Titanium Beta as a Selective Catalyst in the Epoxidation of Bulky Alkenes. Applied Catalysis A: General 1998, 167, 331-342. (4) Morra, E.; Berlier, G.; Borfecchia, E.; Bordiga, S.; Beato, P.; Chiesa, M., Electronic and Geometrical Structure of Zn+ Ions Stabilized in the Porous Structure of Zn-Loaded Zeolite HZsm-5: A Multifrequency Cw and Pulse Epr Study. J. Phys. Chem. C 2017, 121, 14238-14245. (5) Lopez, A.; Soulard, M.; Guth, J. L., Temperature-Induced Monoclinic/Orthorhombic Transition in Germanium Mfi-Type Zeolites. Zeolites 1990, 10, 134-136. (6) Cronstedt, A. F., Natural Zeolite and Minerals. Svenska Vetenskaps Akademiens Handlingar Stockholm 1756, 17, 120. (7)

International Zeolite Association. http://www.iza-online.org/ (accessed 20-6-2017).

(8) Stöcker, M.; Karge, H. G.; Jansen, J. C.; Weitkamp, J., Advanced Zeolite Science and Applications; Elsevier Science, 1994. (9) Davis, M. E., Zeolites and Molecular Sieves: Not Just Ordinary Catalysts. Industrial & Engineering Chemistry Research 1991, 30, 1675-1683. (10) Cacho-Bailo, F.; Etxeberría-Benavides, M.; David, O.; Téllez, C.; Coronas, J., Structural Contraction of Zeolitic Imidazolate Frameworks: Membrane Application on Porous Metallic Hollow Fibers for Gas Separation. ACS Applied Materials & Interfaces 2017, 9, 20787-20796. (11) Wingenfelder, U.; Hansen, C.; Furrer, G.; Schulin, R., Removal of Heavy Metals from Mine Waters by Natural Zeolites. Environmental Science & Technology 2005, 39, 4606-4613. (12) Yeritsyan, H.; Sahakyan, A.; Harutyunyan, V.; Nikoghosyan, S.; Hakhverdyan, E.; Grigoryan, N.; Hovhannisyan, A.; Atoyan, V.; Keheyan, Y.; Rhodes, C., Radiation-Modified Natural Zeolites for Cleaning Liquid Nuclear Waste (Irradiation against Radioactivity). Scientific Reports 2013, 3, 2900. (13) Dyer, A.; Mikhail, K. Y., The Use of Zeolites for the Treatment of Radioactive Waste. Mineralogical Magazine 1985, 49, 203-210. (14) Bish, D. L., Natural Zeolites and Nuclear-Waste Management: The Case of Yucca Mountain, Nevada, USA. In Natural Microporous Materials in Environmental Technology, Misaelides, P.; Macášek, F.; Pinnavaia, T. J.; Colella, C., Eds. Springer Netherlands: Dordrecht, 1999; pp 177-191. 29 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 30 of 38

(15) Farag, H. A. A.; Ezzat, M. M.; Amer, H.; Nashed, A. W., Natural Gas Dehydration by Desiccant Materials. Alexandria Engineering Journal 2011, 50, 431-439. (16) Shirazian, S.; Ashrafizadeh, S. N., Synthesis of Substrate-Modified Lta Zeolite Membranes for Dehydration of Natural Gas. Fuel 2015, 148, 112-119. (17)

Milton, R. M. Drying of Natural Gas by Adsorption. US 3024867 A, 1962.

(18) Wang, S.; Peng, Y., Natural Zeolites as Effective Adsorbents in Water and Wastewater Treatment. Chem. Eng. J. 2010, 156, 11-24. (19) Rožić, M.; Cerjan-Stefanović, Š.; Kurajica, S.; Vančina, V.; Hodžić, E., Ammoniacal Nitrogen Removal from Water by Treatment with Clays and Zeolites. Water Res. 2000, 34, 3675-3681. (20) Tiwari, D. K.; Behari, J.; Sen, P., Application of Nanoparticles in Waste Water Treatment. World Applied Science Journal 2008, 3, 417-433. (21) Gates, B. C., Fluid Catalytic Cracking with Zeolite Catalysts, Paul B. Venuto and E. Thomas Habib, Jr., Marcel Dekker; American Institute of Chemical Engineers: New York, 1980; Vol. 26, p 156. (22) Rahimi, N.; Karimzadeh, R., Catalytic Cracking of Hydrocarbons over Modified Zsm-5 Zeolites to Produce Light Olefins: A Review. Applied Catalysis A: General 2011, 398, 1-17. (23) Songip, A. R.; Masuda, T.; Kuwahara, H.; Hashimoto, K., Test to Screen Catalysts for Reforming Heavy Oil from Waste Plastics. Applied Catalysis B: Environmental 1993, 2, 153164. (24)

Chu, Y. F. Catalytic Reforming with Improved Zeolite Catalysts. US 4927525 A, 1990.

(25) Lane, G. S.; Modica, F. S.; Miller, J. T., Platinum/Zeolite Catalyst for Reforming NHexane: Kinetic and Mechanistic Considerations. J. Catal. 1991, 129, 145-158. (26) Shirazi, L.; Jamshidi, E.; Ghasemi, M. R., The Effect of Si/Al Ratio of Zsm-5 Zeolite on Its Morphology, Acidity and Crystal Size. Cryst. Res. Technol. 2008, 43, 1300-1306. (27) Marinkovic, B. A.; Jardim, P. M.; Rizzo, F.; Saavedra, A.; Lau, L. Y.; Suard, E., Complex Thermal Expansion Properties of Al-Containing Hzsm-5 Zeolite: A X-Ray Diffraction, Neutron Diffraction and Thermogravimetry Study. Microporous Mesoporous Mater. 2008, 111, 110-116. (28) Marinkovic, B. A.; Jardim, P. M.; Saavedra, A.; Lau, L. Y.; Baehtz, C.; de Avillez, R. R.; Rizzo, F., Negative Thermal Expansion in Hydrated Hzsm-5 Orthorhombic Zeolite. Microporous Mesoporous Mater. 2004, 71, 117-124.

30 ACS Paragon Plus Environment

Page 31 of 38 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

(29) Flanigen, E. M.; Bennett, J. M.; Grose, R. W.; Cohen, J. P.; Patton, R. L.; Kirchner, R. M.; Smith, J. V., Silicalite, a New Hydrophobic Crystalline Silica Molecular Sieve. Nature 1978, 271, 512-516. (30) Olson, D. H.; Kokotailo, G. T.; Lawton, S. L.; Meier, W. M., Crystal Structure and Structure-Related Properties of Zsm-5. The Journal of Physical Chemistry 1981, 85, 2238-2243. (31) van Koningsveld, H.; Jansen, J. C.; van Bekkum, H., The Monoclinic Framework Structure of Zeolite H-Zsm-5. Comparison with the Orthorhombic Framework of as-Synthesized Zsm-5. Zeolites 1990, 10, 235-242. (32) van Koningsveld, H.; Tuinstra, F.; Van Bekkum, H.; Jansen, J. C., The Location of PXylene in a Single Crystal of Zeolite H-Zsm-5 with a New, Sorbate‐Induced, Orthorhombic Framework Symmetry. Acta Crystallographica Section B 1989, 45, 423-431. (33) Fyfe, C. A.; Kennedy, G. J.; De Schutter, C. T.; Kokotailo, G. T., Sorbate-Induced Structural Changes in Zsm-5 (Silicalite). J. Chem. Soc., Chem. Commun. 1984, 541-542. (34) Kokotailo, G. T.; Fyfe, C. A.; Kennedy, G. J.; Gobbi, G. C.; Strobl, H.; Pasztor, C. T.; Barlow, G. E.; Bradley, S.; Murphy, W. J.; Ozubko, R. S., Zeolite Structural Investigations by High Resolution Solid State Mas Nmr (Magic Angle Spinning Nuclear Magnetic Resonance). Pure Appl. Chem. 1986, 58, 1367-1374. (35) van Koningsveld, H.; Jansen, J. C.; van Bekkum, H., The Orthorhombic/Monoclinic Transition in Single Crystals of Zeolite Zsm-5. Zeolites 1987, 7, 564-568. (36) Klinowski, J.; Carpenter, T. A.; Gladden, L. F., High-Resolution Solid-State N.M.R. Studies of Temperature-Induced Phase Transitions in Silicalite (Zeolite Zsm—5). Zeolites 1987, 7, 73-78. (37) van Koningsveld, H., High-Temperature (350 K) Orthorhombic Framework Structure of Zeolite H-Zsm-5. Acta Crystallographica Section B 1990, 46, 731-735. (38) Wu, E. L.; Lawton, S. L.; Olson, D. H.; Rohrman, A. C.; Kokotailo, G. T., Zsm-5-Type Materials. Factors Affecting Crystal Symmetry. The Journal of Physical Chemistry 1979, 83, 2777-2781. (39) Quartieri, S.; Arletti, R.; Vezzalini, G.; Di Renzo, F.; Dmitriev, V., Elastic Behavior of Mfi-Type Zeolites: 3 – Compressibility of Silicalite and Mutinaite. J. Solid State Chem. 2012, 191, 201-212. (40) Marra, G. L.; Tozzola, G.; Leofanti, G.; Padovan, M.; Petrini, G.; Genoni, F.; Venturelli, B.; Zecchina, A.; Bordiga, S.; Ricchiardi, G., Orthorhombic and Monoclinic Silicalites: Structure, Morphology, Vibrational Properties and Crystal Defects. In Stud. Surf. Sci. Catal., J. Weitkamp, H. G. K. H. P.; Hölderich, W., Eds. Elsevier: 1994; Vol. Volume 84, pp 559-566.

31 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 32 of 38

(41) Hay, D. G.; Jaeger, H.; West, G. W., Examination of the Monoclinic/Orthorhombic Transition in Silicalite Using Xrd and Silicon Nmr. The Journal of Physical Chemistry 1985, 89, 1070-1072. (42) Hay, D. G.; Jaeger, H., Orthorhombic-Monoclinic Phase Changes in Zsm-5 Zeolite/Silicalite. J. Chem. Soc., Chem. Commun. 1984, 1433-1433. (43) de Vos Burchart, E.; Van Bekkum, H.; Van de Graaf, B., Molecular Mechanics Studies on Mfi-Typezeolites: Part 3. The Monoclinicorthorhombic Phase Transition. Zeolites 1993, 13, 212-215. (44) Haines, J.; Cambon, O.; Levelut, C.; Santoro, M.; Gorelli, F.; Garbarino, G., Deactivation of Pressure-Induced Amorphization in Silicalite Sio2 by Insertion of Guest Species. J. Am. Chem. Soc. 2010, 132, 8860-8861. (45) Sartbaeva, A.; Haines, J.; Cambon, O.; Santoro, M.; Gorelli, F.; Levelut, C.; Garbarino, G.; Wells, S. A., Flexibility Windows and Compression of Monoclinic and Orthorhombic Silicalites. Physical Review B 2012, 85, 064109. (46) Ehrenfest, P., Phasenumwandlungen Im Ueblichen Und Erweiterten Sinn, Classifiziert Nach Den Entsprechenden Singularitaeten Des Thermodynamischen Potentiales; N. V. NoordHollandsche Uitgevers Maatschappij, 1933. (47) Jaeger, G., The Ehrenfest Classification of Phase Transitions: Introduction and Evolution. Archive for History of Exact Sciences 1998, 53, 51-81. (48) Ardit, M.; Martucci, A.; Cruciani, G., Monoclinic–Orthorhombic Phase Transition in Zsm-5 Zeolite: Spontaneous Strain Variation and Thermodynamic Properties. The Journal of Physical Chemistry C 2015, 119, 7351-7359. (49) Mowat, J. P. S.; Miller, S. R.; Griffin, J. M.; Seymour, V. R.; Ashbrook, S. E.; Thompson, S. P.; Fairen-Jimenez, D.; Banu, A.-M.; Düren, T.; Wright, P. A., Structural Chemistry, Monoclinic-to-Orthorhombic Phase Transition, and Co2 Adsorption Behavior of the Small Pore Scandium Terephthalate, Sc2(O2cc6h4co2)3, and Its Nitro- and AminoFunctionalized Derivatives. Inorg. Chem. 2011, 50, 10844-10858. (50) Sharma, S. M.; Sikka, S. K., Pressure Induced Amorphization of Materials. Prog. Mater Sci. 1996, 40, 1-77. (51) Greaves, G. N.; Meneau, F.; Sapelkin, A.; Colyer, L. M.; ap Gwynn, I.; Wade, S.; Sankar, G., The Rheology of Collapsing Zeolites Amorphized by Temperature and Pressure. Nat Mater 2003, 2, 622-629. (52) Neville, G.; Florian, M., Probing the Dynamics of Instability in Zeolitic Materials. J. Phys.: Condens. Matter 2004, 16, S3459. (53) Greaves, G. N.; Meneau, F.; Kargl, F.; Ward, D.; Holliman, P.; Albergamo, F., Zeolite Collapse and Polyamorphism. J. Phys.: Condens. Matter 2007, 19, 415102. 32 ACS Paragon Plus Environment

Page 33 of 38 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

(54) Greaves, G. N.; Meneau, F.; Majérus, O.; Jones, D. G.; Taylor, J., Identifying Vibrations That Destabilize Crystals and Characterize the Glassy State. Science 2005, 308, 1299. (55) Wilding, M. C.; Wilson, M.; McMillan, P. F., Structural Studies and Polymorphism in Amorphous Solids and Liquids at High Pressure. Chem. Soc. Rev. 2006, 35, 964-986. (56) Fu, Y.; Song, Y.; Huang, Y., An Investigation of the Behavior of Completely Siliceous Zeolite Zsm-5 under High External Pressures. The Journal of Physical Chemistry C 2012, 116, 2080-2089. (57) Sartbaeva, A.; Wells, S. A.; Treacy, M. M. J.; Thorpe, M. F., The Flexibility Window in Zeolites. Nat Mater 2006, 5, 962-965. (58) Park, S. H.; Große Kunstleve, R. W.; Graetsch, H.; Gies, H., The Thermal Expansion of the Zeolites Mfi, Afi, Doh, Ddr, and Mtn in Their Calcined and as Synthesized Forms. In Stud. Surf. Sci. Catal., Chon, H.; Ihm, S.-K.; Uh, Y. S., Eds. Elsevier: 1997; Vol. 105, pp 1989-1994. (59) Lassinantti Gualtieri, M.; Gualtieri, A. F.; Hedlund, J., The Influence of Heating Rate on Template Removal in Silicalite-1: An in Situ Ht-Xrpd Study. Microporous Mesoporous Mater. 2006, 89, 1-8. (60) Sen, S.; Wusirika, R. R.; Youngman, R. E., High Temperature Thermal Expansion Behavior of H[Al]Zsm-5 Zeolites: The Role of Brønsted Sites. Microporous Mesoporous Mater. 2006, 87, 217-223. (61) Choi, J.; Jeong, H.-K.; Snyder, M. A.; Stoeger, J. A.; Masel, R. I.; Tsapatsis, M., Grain Boundary Defect Elimination in a Zeolite Membrane by Rapid Thermal Processing. Science 2009, 325, 590. (62) Krokidas, P. G.; Nikolakis, V.; Burganos, V. N., Heating and Sorption Effects on Silicalite-1 Unit Cell Size and Geometry. Microporous Mesoporous Mater. 2012, 155, 65-70. (63) Reisner, B. A.; Lee, Y.; Hanson, J. C.; Jones, G. A.; Parise, J. B.; Corbin, D. R.; Toby, B. H.; Freitag, A.; Larese, J. Z.; Kahlenberg, V., Understanding Negative Thermal Expansion and 'Trap Door' Cation Relocations in Zeolite Rho. Chem. Commun. 2000, 2221-2222. (64) Giddy, A. P.; Dove, M. T.; Pawley, G. S.; Heine, V., The Determination of Rigid-Unit Modes as Potential Soft Modes for Displacive Phase Transitions in Framework Crystal Structures. Acta Crystallographica Section A 1993, 49, 697-703. (65) Bieniok, A.; Hammonds, K. D., Rigid Unit Modes and the Phase Transition and Structural Distortions of Zeolite Rho. Microporous Mesoporous Mater. 1998, 25, 193-200. (66) Evans, J. S. O., Negative Thermal Expansion Materials. J. Chem. Soc., Dalton Trans. 1999, 3317-3326. (67) Heine, V.; Welche, P. R. L.; Dove, M. T., Geometrical Origin and Theory of Negative Thermal Expansion in Framework Structures. J. Am. Ceram. Soc. 1999, 82, 1793-1802. 33 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

(68)

Page 34 of 38

Scienomics Maps. http://www.scienomics.com/ (accessed 20-9-2017).

(69) June, R. L.; Bell, A. T.; Theodorou, D. N., Prediction of Low Occupancy Sorption of Alkanes in Silicalite. The Journal of Physical Chemistry 1990, 94, 1508-1516. (70) June, R. L.; Bell, A. T.; Theodorou, D. N., Molecular Dynamics Study of Methane and Xenon in Silicalite. The Journal of Physical Chemistry 1990, 94, 8232-8240. (71) June, R. L.; Bell, A. T.; Theodorou, D. N., Transition-State Studies of Xenon and Sulfur Hexafluoride Diffusion in Silicalite. The Journal of Physical Chemistry 1991, 95, 8866-8878. (72) Goodbody, S. J.; Watanabe, K.; MacGowan, D.; Walton, J. P. R. B.; Quirke, N., Molecular Simulation of Methane and Butane in Silicalite. J. Chem. Soc., Faraday Trans. 1991, 87, 1951-1958. (73) Kiselev, A. V.; Lopatkin, A. A.; Shulga, A. A., Molecular Statistical Calculation of Gas Adsorption by Silicalite. Zeolites 1985, 5, 261-267. (74) June, R. L.; Bell, A. T.; Theodorou, D. N., Molecular Dynamics Studies of Butane and Hexane in Silicalite. The Journal of Physical Chemistry 1992, 96, 1051-1060. (75) Maginn, E. J.; Bell, A. T.; Theodorou, D. N., Transport Diffusivity of Methane in Silicalite from Equilibrium and Nonequilibrium Simulations. The Journal of Physical Chemistry 1993, 97, 4173-4181. (76) Snurr, R. Q.; June, R. L.; Bell, A. T.; Theodorou, D. N., Molecular Simulations of Methane Adsorption in Silicalite. Molecular Simulation 1991, 8, 73-92. (77) Makrodimitris, K.; Papadopoulos, G. K.; Theodorou, D. N., Prediction of Permeation Properties of Co2 and N2 through Silicalite Via Molecular Simulations. The Journal of Physical Chemistry B 2001, 105, 777-788. (78) O’Malley, A. J.; Catlow, C. R. A.; Monkenbusch, M.; Jobic, H., Diffusion of Isobutane in Silicalite: A Neutron Spin–Echo and Molecular Dynamics Simulation Study. The Journal of Physical Chemistry C 2015, 119, 26999-27006. (79) O'Malley, A. J.; Catlow, C. R. A., Molecular Dynamics Simulations of Longer NAlkanes in Silicalite: A Comparison of Framework and Hydrocarbon Models. PCCP 2013, 15, 19024-19030. (80) Kolokathis, P. D.; Pantatosaki, E.; Gatsiou, C.-A.; Jobic, H.; Papadopoulos, G. K.; Theodorou, D. N., Dimensionality Reduction of Free Energy Profiles of Benzene in Silicalite-1: Calculation of Diffusion Coefficients Using Transition State Theory. Molecular Simulation 2014, 40, 80-100. (81) Skoulidas, A. I.; Sholl, D. S., Molecular Dynamics Simulations of Self-Diffusivities, Corrected Diffusivities, and Transport Diffusivities of Light Gases in Four Silica Zeolites to

34 ACS Paragon Plus Environment

Page 35 of 38 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

Assess Influences of Pore Shape and Connectivity. The Journal of Physical Chemistry A 2003, 107, 10132-10141. (82) Krokidas, P.; Castier, M.; Economou, I. G., Computational Study of Zif-8 and Zif-67 Performance for Separation of Gas Mixtures. The Journal of Physical Chemistry C 2017, 121, 17999-18011. (83) Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1-19. (84) Parrinello, M.; Rahman, A., Crystal Structure and Pair Potentials: A Molecular-Dynamics Study. Phys. Rev. Lett. 1980, 45, 1196-1199. (85) Parrinello, M.; Rahman, A., Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182-7190. (86)

The Amber Force Field. http://ambermd.org/#ff (accessed 24-10-17).

(87) Mayo, S. L.; Olafson, B. D.; Goddard, W. A., Dreiding: A Generic Force Field for Molecular Simulations. The Journal of Physical Chemistry 1990, 94, 8897-8909. (88) Sun, H., Compass:  An Ab Initio Force-Field Optimized for Condensed-Phase Applicationsoverview with Details on Alkane and Benzene Compounds. The Journal of Physical Chemistry B 1998, 102, 7338-7364. (89) Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T., Structure and Energetics of Ligand Binding to Proteins: Escherichia Coli Dihydrofolate Reductase-Trimethoprim, a Drug-Receptor System. Proteins: Structure, Function, and Bioinformatics 1988, 4, 31-47. (90) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M., Uff, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024-10035. (91) Hill, J. R.; Sauer, J., Molecular Mechanics Potential for Silica and Zeolite Catalysts Based on Ab Initio Calculations. 1. Dense and Microporous Silica. The Journal of Physical Chemistry 1994, 98, 1238-1244. (92) Boulfelfel, S. E.; Ravikovitch, P. I.; Koziol, L.; Sholl, D. S., Improved Hill–Sauer Force Field for Accurate Description of Pores in 8-Ring Zeolites. The Journal of Physical Chemistry C 2016, 120, 14140-14148. (93) Nosé, S., A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. The Journal of Chemical Physics 1984, 81, 511-519. (94) Hoover, W. G., Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695-1697.

35 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 36 of 38

(95) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R., A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. The Journal of Chemical Physics 1982, 76, 637-649. (96) Tsalikis, D. G.; Koukoulas, T.; Mavrantzas, V. G.; Pasquino, R.; Vlassopoulos, D.; Pyckhout-Hintzen, W.; Wischnewski, A.; Monkenbusch, M.; Richter, D., Microscopic Structure, Conformation, and Dynamics of Ring and Linear Poly(Ethylene Oxide) Melts from Detailed Atomistic Molecular Dynamics Simulations: Dependence on Chain Length and Direct Comparison with Experimental Data. Macromolecules 2017, 50, 2565-2584. (97) Leardini, L.; Quartieri, S.; Vezzalini, G., Compressibility of Microporous Materials with Cha Topology: 1. Natural Chabazite and Sapo-34. Microporous Mesoporous Mater. 2010, 127, 219-227. (98) Lindemann, F. A., Über Die Berechnung Molekularer Eigenfrequenzen. Phys. Z. 1910, 11, 609-612. (99) Kaiwang, Z.; Stocks, G. M.; Jianxin, Z., Melting and Premelting of Carbon Nanotubes. Nanotechnology 2007, 18, 285703. (100) Koparde, V. N.; Cummings, P. T., Sintering of Titanium Dioxide Nanoparticles: A Comparison between Molecular Dynamics and Phenomenological Modeling. J. Nanopart. Res. 2008, 10, 1169-1182. (101) Lempesis, N.; in ‘t Veld, P. J.; Rutledge, G. C., Simulation of the Structure and Mechanics of Crystalline 4,4′-Diphenylmethane Diisocyanate (Mdi) with N-Butanediol (Bdo) as Chain Extender. Polymer 2016, 107, 233-239. (102) Goudeli, E.; Pratsinis, S. E., Surface Composition and Crystallinity of Coalescing Silver– Gold Nanoparticles. ACS Nano 2017, 11, 11653-11660. (103) Yang, X.; Gan, C.; Xiong, H.; Huang, L.; Luo, X., Fabrication and Characterization of Sio2@Tio2@Silicalite-1 Catalyst and Its Application for Degradation of Rhodamine B. RSC Advances 2016, 6, 105737-105743. (104) Kestell, J. D.; Zhong, J.-Q.; Shete, M.; Waluyo, I.; Sadowski, J. T.; Stacchiola, D. J.; Tsapatsis, M.; Boscoboinik, J. A., Studying Two-Dimensional Zeolites with the Tools of Surface Science: Mfi Nanosheets on Au(111). Catal. Today 2017, 280, Part 2, 283-288. (105) Elyassi, B.; Jeon, M. Y.; Tsapatsis, M.; Narasimharao, K.; Basahel, S. N.; Al-Thabaiti, S., Ethanol/Water Mixture Pervaporation Performance of B-Oriented Silicalite-1 Membranes Made by Gel-Free Secondary Growth. AlChE J. 2016, 62, 556-563. (106) Lovallo, M. C.; Gouzinis, A.; Tsapatsis, M., Synthesis and Characterization of Oriented Mfi Membranes Prepared by Secondary Growth. AlChE J. 1998, 44, 1903-1913.

36 ACS Paragon Plus Environment

Page 37 of 38 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

(107) Torres, C.; Gulín-González, J.; Navas-Conyedo, E.; Demontis, P.; Suffritti, G. B., The Behavior of Silicalite-1 under High Pressure Conditions Studied by Computational Simulation. Struct. Chem. 2013, 24, 909-915. (108) Waheed, N.; Lavine, M. S.; Rutledge, G. C., Molecular Simulation of Crystal Growth in N-Eicosane. The Journal of Chemical Physics 2002, 116, 2301-2309.

37 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 38 of 38

TOC graphic

38 ACS Paragon Plus Environment