Improved HillâSauer Force Field for Accurate ... - ACS Publications

Subscriber access provided by UNIV OF NEBRASKA - LINCOLN

Article

An Improved Hill-Sauer Force Field for Accurate Description of Pores in 8-Membered Ring Zeolites Salah Eddine Boulfelfel, Peter I. Ravikovitch, Lucas Koziol, and David S. Sholl J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b03674 • Publication Date (Web): 10 Jun 2016 Downloaded from http://pubs.acs.org on June 15, 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 21

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

An Improved Hill-Sauer Force Field for Accurate Description of Pores in 8-Ring Zeolites Salah Eddine Boulfelfela, Peter I. Ravikovitchb*, Lucas Koziolb, and David S. Sholla** a

School of Chemical and Biomolecular Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0100, United States b Corporate Strategic Research, ExxonMobil Research and Engineering, 1545 Route 22 East, Annandale, New Jersey 08801, United States * +1-908-730-2280 [email protected] ** +1-404-894-2822 [email protected]

Abstract We have introduced a simple modification of the well-known Hill-Sauer force field for silica. The modified force field improves the accuracy with which pore sizes and framework flexibility in small pore zeolites are described. The modification focused on the Si–O–Si and O–Si–O angles in these materials, which are instrumental in controlling vibrations of the framework such as twisting of the near-rigid SiO4 units. The accuracy of the modified HillSauer force field was compared with data from extensive density functional theory calculations of zeolite structures and dynamics. The transferability of the force field was tested on 13 experimentally known silica 8-ring zeolites. Additional tests examining the thermal expansion behavior of selected zeolites were also performed.

1 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 21

1. Introduction Eight-membered ring (8-ring) zeolites such as DDR, CHA, LTA, and ERI are an important class of small-pore zeolites with effective pore sizes between 3.4–4.6 Å1-2. They are promising candidates for industrial processes including ion exchange, heterogeneous catalysis, and mixture separations3-5. Figure 1 illustrates the 8-ring in the silica form of LTA. The window size in 8-ring zeolites is comparable to the kinetic diameter of small hydrocarbons, which makes them potential candidates for the separation of mixtures like ethane–ethene, propane–propylene, and isomers of butane1, 6-8. The small differences in the pore dimensions among various 8-ring zeolites enable the molecular sieving properties to be adjusted in a controllable fashion1. For example, Hedin et al. investigated the effect of window size on diffusion in LTA structures in a series of experiments9. They studied the diffusion of propylene in CaA, ITQ-29, and partially exchanged NaCaA. The zeolites LTA5A and ITQ-29 have smaller 8-ring windows unobstructed by cations but showed diffusivity lower by an order of magnitude compared to the cationic structures CaA and partially exchanged NaCaA. The efficiency of hydrocarbon separation processes in 8-ring zeolites relies on a combination of adsorption and diffusion selectivity of the molecules in the target mixture2, 10. In addition to the important role of pore size in the sieving properties of 8-ring zeolites, the flexibility of the framework, especially the atoms forming the window, is another factor controlling molecular diffusion in these materials. A number of modeling studies have shown that characteristics of the 8-ring window due to intra-framework vibrations and related motions are critical to quantitatively understanding molecular diffusion in these materials11-12. In terms of theoretical simulations, the ability to correctly predict the pore size and framework flexibility of 8-ring zeolites relies on the development of reliable force fields 13. The widespread use of siliceous materials either in their dense forms or as porous media has prompted the need for simulations to predict and understand different physical and chemical properties and characteristics of these materials13-16. The literature features numerous force fields for siliceous materials13, 16-28. Hill and Sauer17-18 developed one of the most widely used force fields for the simulation of pure silica zeolites. It is often used to compute properties of zeolites ranging from the simulation of diffusion process with flexible frameworks2, 11-12, 29-30 to spectroscopic and vibrational studies of zeolite materials14, 31. Figure 2 shows the window size distribution in pure silica LTA computed using the HillSauer force field and also from first principles molecular dynamics at the DFT level using the PBE functional. Further details on the methods used for these calculations are given below. The window size observed experimentally32 is also shown in Figure 2. Although the bimodal nature of the window size distribution is reproduced with both approaches, the mean window sizes for WS1 (WS2) from the Hill-Sauer force field are shifted in comparison to experiments by ~−0.4 Å (~+0.4 Å). Differences of this magnitude can have a large impact on the calculations of diffusivity when the adsorbate size is comparable to that of the window1, 12. The results in Figure 2 strongly suggest that the Hill-Sauer force field is not ideally suited to accurately describing framework flexibility in 8-ring zeolites.


Page 3 of 21

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. A ball-and-stick representation of pure silica LTA unit cell with oxygen and silicon ions in red and blue, respectively. Oxygen ions numbered from 1 to 8 define the pore with two characteristic distances WS1 (1–5 black stick) and WS2 (2–6 red stick).

Figure 2. Window size distribution in LTA sampled from NpT molecular dynamics simulations using Hill-Sauer force field (HSFF) and DFT (PBE functional). In each case, the distributions for WS1 and WS2 (defined in Figure 1) are shown separately. Solid and dashed blue lines correspond to experimental values of WS1 and WS2, respectively. 3 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 4 of 21

In this paper, we introduce a modified version of the Hill-Sauer force field (HSFF) that allows an improved description of pore size and framework flexibility in 8-ring zeolites. The modifications focus on the angle terms while keeping the bonding terms of the rigid SiO4 tetrahedral units intact. We tested the overall transferability of the modified version of the HSFF across a series of experimentally-accessible 8-ring zeolites. We find that in addition to better describing 8-ring windows, the prediction of thermal expansion behavior of 8-ring zeolites by our modified force field is also improved relative to the original force field. The modified force field we introduce below should be useful in efforts to quantitatively predict the properties of small pore zeolites for molecular separations and other applications.

2. Simulation Details 2.1. First-principles Molecular Dynamics and Structural Optimization To provide reliable information on the window size distribution in 8-ring zeolites, BornOppenheimer molecular dynamics was carried out using CP2K (version 2.6)33. Energies and forces were computed from density functional theory (DFT) as implemented in the module QUICKSTEP34. In these calculations, the self-consistent field (SCF) minimizer was based on the orbital transformation method35, and a mixed Gaussian and Plane-Wave (GPW) method3637 was used in combination with Goedecker-Teter-Hutter (GTH) pseudopotentials38-39. The plane wave and DZVP-MOLOPT-SR-GTH auxiliary basis sets cutoff were 850 and 70 Ry, respectively. We report below on calculations using the local density approximation (LDA)38 and two gradient-corrected functionals, BLYP40-42, and PBE43. First-principles Molecular Dynamics (MD) simulations were propagated at 300 K and 1 bar in the NpT ensemble using the Nosé-Hoover thermostat44 and a time-step of 1 fs. Prior to MD, full-cell structural optimization of 8-ring zeolites was carried out using the same settings as for first-principles molecular dynamics using DZVP-MOLOPT-GTH instead of DZVP-MOLOPT-SR-GTH as the auxiliary basis set. Geometry and cell optimization were performed simultaneously with the stress tensor computed at each step using a conjugate gradient methods (CG) for the cell and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method for geometry. Convergence criterion for maximum geometry and root-mean-square geometry changes were set to 3.0 10-3 and 1.5 10-3 Bohr. Convergence criterion for maximum force component and root-meansquare force changes were set to 3.0 10-4 and 1.5 10-4 Bohr-1Hartree. In both molecular dynamics and structural optimization we used a 2 × 2 × 2 super-cell of LTA and CHA (trigonal setting).

2.2 Classical Molecular Dynamics and Structural Optimization Classical molecular dynamics was carried out using the LAMMPS code45. Simulations were performed at 300 K and 1 bar in the NpT ensemble using a Nosé-Hoover thermostat44 and a time-step set of 1 fs. The Hill-Sauer force field17-18 was used to describe van der Waals, bond, angle, dihedral, torsion, and cross interactions. Long range electrostatic interactions were computed using the Ewald method44. For accurate force field fitting a cutoff equal to 15 Å was used for electrostatics and van der Waals interactions. For production simulations a smaller cutoff for van der Waals interactions can be used. Prior to MD, a conjugate gradient minimizer was employed for structural optimization including cell shape and volume and atomic positions. Stopping tolerance criteria for energy and forces were 10-8 and 10-10 kcal.mol-1/Å, respectively. In both molecular dynamics and structural optimization we used a 2 × 2 × 2 super-cell of LTA and CHA (hexagonal setting). 4 ACS Paragon Plus Environment

Page 5 of 21

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


3. Results 3.1. Force-field Fitting The Hill-Sauer force field can be written as:

VHS =Vbond +Vangle +Vtorsion+Vbond-bond +Vangle-angle +Vbond-angle +Vbond-angle-torsion+VLJ +VCoulomb 1 A full description of different terms and corresponding parameters is reported in section S1 of the Supporting Information. Pure silica zeolites form varied topologies based on different arrangements of SiO4 building units. SiO4 units are known to be relatively rigid and the flexibility of zeolite frameworks is largely due to vibrational properties of Si–O–Si and O–Si–O angles. The angle terms in the HSFF are expressed as:

Vangle =K 2 − ! 2 +K 3 − ! 3 +K 4 − ! 4 2 In the original HSFF, ! , K $ , K 3 , and K % are 112.02°, 81.9691 kcal/mol/rad2, -36.5814 kcal/mol/rad3, and 116.9558 kcal/mol/rad4 for O–Si–O and 173.7651°, 20.7015 kcal/mol/rad2, 27.5506 kcal/mol/rad3, and 10.9930 kcal/mol/rad4 for Si–O–Si. We modified the HSFF by tuning the equilibrium angle ! for Si–O–Si and O–Si–O while keeping the constants K $ , K & , and K % unchanged. These constants were obtained from ab initio calculations with an accurate basis set: a double zeta with polarization (DZP) and a triple zeta with polarization (TZP) for Si and O atoms, respectively. However, equilibrium angles are not adequate for an accurate description of bulk pure silica zeolites because they were determined from calculations on small Si-O models (molecular systems and clusters). The tuning of the equilibrium angles, ! , was done in two steps by focusing on the LTA structure as a prototypical 8-ring zeolite. First we fixed ! for the O–Si–O angle and performed a series of full cell optimizations at different values of ! for the Si–O–Si angle (Figure 3, top panel). The Si–O–Si angle corresponding to the best agreement between calculated and experimental32 WS1 (3.995 Å) and WS2 (4.218 Å) was selected as the new ! (Si–O–Si), giving a value of 150°. Subsequently ! for the Si–O–Si angle was fixed to this value while another series of full cell optimizations was performed as a function of ! for the O–Si–O angle (Figure 3, lower panel). As shown in Figure 3, setting the O–Si–O angle to 113° gives the best agreement between the calculated and experimental WS1 and WS2. As a result, this choice was set as new equilibrium ! for this angle. Figure 4 shows the O–Si–O and Si–O–Si angle terms in the original and modified HSFF. There is only a minor difference between the two potentials for the O–Si–O angle, while the modified Si–O–Si potential features a shift to lower values of Si–O–Si angle. By construction, the shape of the two potential is unchanged in the modified HSFF; only the equilibrium angle was changed. Although the strategy behind this modification of the HSFF is simple, we show below that it has a great impact on the calculated pore size and the dynamical behavior of 8-ring windows in zeolites.



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 21

Figure 3. Variation of window sizes WS1 (circles) and WS2 (squares) as a function of Si–O– Si and O–Si–O equilibrium angles in the 8-ring window of LTA. Dashed black and dotted red lines correspond to the experimental values of window sizes WS1 and WS2, respectively. 3.2. Window Size Distributions In this section, a detailed description of the pore size and dynamical behavior of 8-ring windows is presented for ITQ-29 (the siliceous form of LTA) and chabazite (CHA). All calculations are for pure silica materials. A comparison between the original and modified HSFF is shown for both structures and validated using data collected from first-principles molecular dynamics. This analysis is further extended to cover 13 experimentally known 8ring zeolites to test the transferability of the modified HSFF. Pore size or window size distributions, WSD, are defined as T

N

WSD= ) ) dij t 3 t=0 w=1


Page 7 of 21

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


For the WS1 distribution, (i,j)=(1,5) and (3,7), while for the WS2 distribution, (i,j)=(2,6) and (4,8) (see Figure 1 for definition of indices used to label atoms in the 8-ring). Here, t, T, w, and N are simulation time, total simulation time, window index, and total number of windows, respectively. The window size distributions in LTA and CHA were collected by recording two characteristic distances, O1–O5 and O3–O7 for WS1 and O2–O6 and O4–O8 for WS2, for all 8-ring windows in the simulation volume during the molecular dynamics simulation after subtracting one oxygen ion diameter (2.7 Å) from the O–O distances46. The total window size distribution was defined as the sum of WS1 and WS2.

Figure 4. The Si–O–Si and O–Si–O angle contributions to the Hill-Sauer force field using the original (solid lines) and modified (dashed lines) equilibrium angles. In addition to the window size distributions defined above, minimal window size distributions were also computed for LTA and CHA. Experiments and simulations1-2, 11-12 have shown that the minimum free diameter of the window in 8-ring zeolites is a critical dimension for important processes such as diffusion. Here, we define the minimal window size distribution as the ensemble average of the shortest O-O distance for each 8-ring window,



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

T

Page 8 of 21

N

MWSD= ) ) min0d15 t,d37 t,d26 t,d48 t6 4 t=0 w=1

3.2.1 LTA We first consider the window size distributions for LTA. Figure 5 compares the window size distribution obtained using the original and modified HSFF. Both versions of the force field reproduce the bimodal distribution previously reported in theoretical investigations of LTA structures using different force fields2, 30. Although the window size distributions for WS1 and WS2 from the modified HSFF show a significant overlap, there is a clear separation between the distributions when using the original HSFF. As expected because of the approach we took to developing the HSFF, the mean value of the window size distributions from the modified HSFF are in excellent agreement with pore sizes from experiments32. As already shown in Figure 1, the window size distributions computed using the original HSFF underestimate/overestimate WS1/WS2 by ~0.4 Å. For production NVT molecular dynamics simulations it is important to perform the calculations at the equilibrium volume predicted by HSFF and not the experimental value (see Figure S1 in Supporting Information). Experimentally, only the mean values of the window size distributions shown in Figure 5 are available. In order to further examine the quality of the modified HSFF we compared data obtained from classical force fields to that obtained from DFT using several functionals. Figure 5 shows a comparison between window size distributions obtained from firstprinciples NpT molecular dynamics simulations and using the modified HSFF. Results using DFT-LDA and DFT-PBE are shown in Figure 5 and DFT-BLYP in Figure S2 in Supporting Information. Although there are variations in the results among the different functionals, many of the characteristics of the window size distributions are the same in each set of DFT data. There is a significant overlap between the WS1 and WS2 distributions in each set of DFT calculations. In terms of the mean pore size, the DFT-PBE results are in better agreement with experimental data (and the modified HSFF) than DFT-BLYP and DFT-LDA. All of the DFT calculations give somewhat broader window size distributions than the modified HSFF. For example, DFT-PBE gives a standard deviation of 0.196 Å (0.160 Å) for WS1 (WS2), while the modified HSFF gives 0.115 Å (0.099 Å). We also examined the validity of the modified HSFF by assessing the temperature dependence of the window size distributions. Both the total and minimal window size distributions were analyzed at 300, 600, 900, and 1200 K. The comparison between the total window size distributions using the original and modified HSFF is shown in Figure S2 in the Supporting Information. Upon increasing temperature, a widening of the total pore size distribution was observed with an increased overlap between WS1 and WS2 distributions for the original HSFF. The bimodal distribution for the total pore size distribution becomes less distinct as temperature is increased when using the modified HSFF. A comparison between the minimal window size distributions at different temperatures using the original and modified HSFF is shown in Figure 6. While the mean value of this distribution is almost independent of temperature using the original HSFF, a decrease in the mean minimal pore size with temperature is observed with the modified HSFF. To further probe this significant difference in temperature dependence, first principles molecular 8 ACS Paragon Plus Environment

Page 9 of 21

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


dynamics simulations using the same three functionals as above (PBE and LDA in Figure 6 and BLYP in Figure S3 in Supporting Information) were performed. Figure 6 shows there is good agreement between the temperature dependence of the minimal pore size distribution predicted by the modified HSFF and the DFT calculations; the shift of the distribution predicted by the modified HSFF is clearly seen in each of the DFT calculations. As in Figure 5, the distributions observed from DFT are wider than from our classical simulations.

Figure 5. Window size distributions in LTA sampled from classical and first-principles NpT molecular dynamics simulations. In each case, the distributions for WS1 and WS2 (defined in Figure 1) are shown separately. Solid and dashed blue lines indicate the experimentally observed WS1 and WS2, respectively.



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 21

Figure 6. Temperature dependence of minimal window size distributions in LTA sampled from classical and first-principles NpT molecular dynamics simulations. 3.2.2 CHA The calculations presented above for LTA were also performed for the all silica form of CHA. The mean window sizes WS1 and WS2 computed using the original HSFF show an apparent agreement with experimental data47, as shown on Figure 7. A surprising feature of this result is that the distribution for WS1 shows a bimodal distribution. If this distribution is resolved around two mean WS distances rather than one, the two means, WS1a and WS1b, correspond to 3.45 and 4.15 Å, respectively. The bimodal distribution for WS1 was not observed when the modified HSFF was used. The mean pore sizes of WS1 and WS2 (3.80 and 4.25 Å) from the modified HSFF are in a good agreement with experimental values (3.77 and 4.23 Å)2, 47. Unlike the similar result above for LTA, this can be viewed as a test of the modified force field because this information was not included in fitting the potential. The window size distributions in CHA were obtained from first-principles NpT molecular dynamics using three different functionals (PBE and LDA in Figure 7 and BLYP in Figure S4 in Supporting Information). The bimodal distribution predicted by the original HSFF was 10 ACS Paragon Plus Environment

Page 11 of 21

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


not observed in any of our DFT calculations. The window size distributions observed using DFT are similar to the results from the modified HSFF, although as in the case of LTA DFT predicts wider distributions and more significant overlap between the WS1 and WS2 distributions than in the classical simulations. The temperature dependence of minimal window size distribution was also computed (see Figure S5 in the Supporting Information). A decrease in mean minimal window size was observed upon increasing temperature when the modified HSFF was used. The original HSFF showed the inverse trend with wider distributions. Experimentally48-49, the volume of CHA decreased when temperature was increased from 250 to 875 K. Carey et al50-51 directly linked the decrease in volume with a mechanism of expansion in window size in 8-ring zeolites. Therefore, the widening of window size distribution resulting into a decrease in minimal window size observed in CHA using modified HSFF is more consistent with experimental observations48-51.

Figure 7. Window size distributions in CHA sampled from classical and first-principles NpT molecular dynamics simulations. In each case, the distributions for WS1 and WS2 (defined in Figure 1) are shown separately. Solid and dashed blue lines indicate experimental WS1 and WS2, respectively. 11 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 21

This analysis was further extended to cover 13 experimentally known 8-ring zeolites (listed in Table S1 in the Supporting Information) to test the transferability of the modified HSFF. The accuracy of both original and modified versions of Hill-Sauer force field is discussed in comparison to DFT results using three functionals: BLYP, LDA, and PBE. To avoid encumbrance of parity plot figures and improve their clarity, results for DFT-BLYP were omitted in Figure 8, 9, and 10, and reported in SI in Figure S6, S7, and S8. Using fullstructural optimization we computed the two distances characteristic of 8-ring pores, WS1 and WS2, for each material. The comparison to experimental data is shown in Figure 8 and 9 for WS1 and WS2, respectively. The corresponding numerical values are listed in Table S2 in the Supporting Information. Using the original HSFF, the mean absolute error was 0.550 Å (0.436 Å) for WS1 (WS2). The modified HSFF showed a better agreement with experiments as the mean absolute error on WS1 (WS2) was reduced to 0.434 Å (0.252 Å). The performance of the modified HSFF is comparable to the DFT LDA and PBE levels of theory shown in Figure 8 and 9. Results including DFT-BLYP calculations are reported in the Supporting Information in Figure S6 and S7 for WS1 and WS2, respectively. Mean absolute errors on WS1 (WS2) were 0.424 Å (0.248 Å), 0.538 Å (0.240 Å), and 0.439 Å (0.223 Å) using BLYP, LDA, and PBE, respectively.

Figure 8. Window size WS1 in 13 pure silica 8-ring zeolites computed from fully optimized (geometry and cell) structures using HSFF and DFT compared to experimental data.


Page 13 of 21

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


We also determined the energy minimized volume of the same set of 13 8-ring zeolites, as shown in Figure 10 (detailed information on crystallographic cell parameters and angles for each structure is listed in Table S3 in the Supporting Information). Due to large range of volumes, only 9 structures are shown for a better clarity. The mean absolute error was 446.2 Å3 using the original HSFF and 344.8 Å3 using the modified force field. The improved accuracy is comparable to DFT-BLYP calculations with a mean absolute error of 279.8 Å3 (see Figure S8 in the Supporting Information). LDA and PBE based calculations showed smaller mean absolute errors on volume; 66.4 Å3 and 189.2 Å3, respectively. Although deviations between experimentally reported structures and these optimized structures exist, the results in Figure 8-10 indicate that the modified HSFF gives a systematic improvement in the prediction of the structure of 8-ring zeolites relative to the original HSFF.

Figure 9. Window size WS2 in pure silica 8-ring zeolites computed from fully optimized (geometry and cell) structures using HSFF and DFT compared to experiments.



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 21

Figure 10. Volume of pure silica 8-ring zeolites computed from fully optimized (geometry and cell) structures using HSFF and DFT compared to experiments. 3.3. Volume and Thermal Expansion Another use of classical force fields for silica zeolites is to predict the temperature dependence of the unit cell volumes. Capturing the correct change in volume and thermal expansion for 8-ring zeolites can have a direct impact on the quality of prediction of pore size and 8-ring window flexibility. The average volume of LTA and CHA was computed using the original and modified HSFF over a range of temperatures from 100 to 1200 K with a 100 K increment. For each temperature, the simulation system was propagated using the NpT ensemble for 5 ns after 1 ns of equilibration. The corresponding volume was averaged over the production run. As shown in Figure 11, the modified HSFF predicts a smaller volume than the original version for all simulated temperatures and for both structures. Overall, both force fields predict a decrease in volume over the wide temperature range shown in Figure 11. However, original HSFF showed an increase in volume over the range 100–300 K for LTA and 100–200 K for CHA (Figure 11, circle symbols), which is inconsistent with experimental observations48-51. The modified HSFF restored the correct volume change upon temperature increase from 100 to 300 K for LTA and 100 to 200 K for CHA (Figure 11, square symbols). Comparison to experiments of Carey et al.50-51 and Woodcock et al.49 is provided in Figure S9 and S10 in Supporting Information section S7. 14 ACS Paragon Plus Environment

Page 15 of 21

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


Experimentally, temperature dependence of volume of silica zeolites can be evaluated by measuring the thermal expansion coefficient over a range of temperature. This coefficient is defined as: 7v =

VT2 -VT1 5 VT1 T2 -T1

The LTA and CHA structures are known to have the largest negative thermal expansion and volume contraction upon decrease of temperature among silica zeolites48-51. A comparison between thermal expansion coefficient from experiments and simulations using original and modified HSFF is given in Table 1. The negative sign is indicator of volume compression. The inaccurate description of change in volume over the range 100–300 K using original HSFF is reflected by the positive sign of thermal expansion coefficient for LTA and CHA. The coefficients computed using modified HSFF, although different by one order of magnitude from experiments, are predicted with the correct sign over different ranges of temperatures as shown in Table 1. The modification of angle terms in HSFF to improve the description of pore size in 8-ring zeolites was motivated in part by the experimental observations linking Si–O–Si and O–Si–O angles variation and pore size in some zeolites4951 . To include unit cell volume of 8-ring zeolites in the force field fitting procedure, more parameters should be tuned and more terms should be revised (e.g. van der Waals interactions). We note that the discussion above in Figure 10 compared minimum energy structures from various levels of theory (i.e. volumes at a temperature of 0 K) to experimental data that is primarily obtained at room temperature. Figure 11 shows that including the effects of thermal expansion or contraction in that comparison would not change our discussion of Figure 10 in a significant way because the volume changes due to thermal expansion are small relative to the typical differences observed between the calculated and experimental volumes. As a final test of the modified HSFF, we computed the vibrational density of states (VDOS) for LTA using classical molecular dynamics with the original and modified HSFF. The resulting VDOS are shown in on Figure S11 in the Supporting Information. Although there are are of course quantitative differences between the VDOS with the two different FFs, the overall features are similar.



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 21

Figure 11. Temperature dependence of volume in (a) LTA and (b) CHA using original and modified Hill-Sauer force field. Each data point was obtained by averaging system volume over simulation time during NpT MD runs. CHA +0.18 10 (100–300 K) αv (Original Hill-Sauer FF) [K ] −3.37 10-6 (300–900 K) αv (Modified Hill-Sauer FF) [K-1] −2.97 10-6 (100–300 K) −2.75 10-6 (300–900 K) αv (Experiments) [K-1] −26.1 10-6 (293–873 K)50 −16.7 10-6 (293–873 K)49 Table 1. Thermal expansion coefficients for LTA and CHA from theory and experiments. -1

LTA +1.46 10 (100–300 K) −2.21 10-6 (300–900 K) −4.64 10-6 (100–300 K) −3.07 10-6 (300–900 K) −22.1 10-6 (100–300 K)51 -6

-6

4. Conclusions We have introduced a simple modification of the well-known Hill-Sauer force field for silica that was developed to more accurately describe pore sizes and framework flexibility in small pore zeolites. Our modification focused on the Si–O–Si and O–Si–O angles in these materials, which are instrumental in controlling vibrations of the framework such as twisting of the near-rigid SiO4 units. Although we focused on 8-ring zeolites, it is useful to compare the contents of the modified FF with a broader class of silica zeolites. Wragg et al.52 analyzed the structure of 35 well-defined pure-silica zeolites. They reported that the O–Si–O angles ranged from 96.8° to 129.2° with a mean value of 109.5° and a median value of 110°. The Si–O–Si angle, however, adopted a wider range from 133.6° to 180°, with mean value of 154° and median value of 148°. In this work, we used optimized equilibrium values for O– Si–O and Si–O–Si angles of 113° and 150°, respectively, values that are very similar to the median values reported by Wragg et al. This observation supports the relatively good transferability of the modified HSFF we observed for a set of 13 pure-silica 8-ring zeolites.


Page 17 of 21

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


Many of results we have presented focused on LTA and CHA. The change in window size in moving from the original to the modified HSFF for LTA is, at one level, small (~ 0.4 Å). We emphasize, however, that changes in window diameter of this magnitude can have a large impact on the diffusion of molecules through the zeolite windows. Using molecular dynamics simulations, Krishna and van Baten2 showed that a change in pore size of as little as 0.5 Å resulted in a change of about two orders of magnitude in diffusivities of CH4 in 8-ring zeolites at 500 K. It is reasonable to think that the impact of window size will be even more dramatic when the diffusion of larger or more complex molecules is considered11-12, 53-54. Our calculations indicate that the modified HSFF gives an improved description of several subtle features of small pore zeolites compared to the original HSFF. For example, the modified FF better reproduces the temperature-dependent trends in the window size distributions as determined by comparisons with DFT-based calculations using multiple functionals. It seems likely that this property will allow for more precise predictions of the temperature-dependence of molecular diffusion rates in these materials than could be made with the original HSFF. In much of this paper we have considered the accuracy of classical FF calculations via comparisons with data from DFT calculations. It is therefore useful to explicitly note that important problems remain to which DFT-based calculations cannot be practically applied and for which classical force fields are well suited. Molecular diffusion provides an example of this kind of problem. Although approximate methods exist that neglect coupling between diffusing molecules and flexible frameworks11, 54-55, quantitatively predicting the diffusion rates of molecules that are complex or large relative to pore openings intrinsically involves an interplay between molecules near pore mouths and deformations of the pores induced by these molecules12, 53. Calculations based on classical force fields are well suited to analyze this situation. Using DFT-based methods for the same problem, is problematic both because of the very large number of configurations that must be examined and also because of the challenges of accurately describing energetic interactions of physisorbed molecules with pore frameworks using DFT. 5. Supporting Information Hill-Sauer force field description. List and structural details of 13 pure silica 8-ring zeolites. Additional window size distribution and volume data. Experimental and calculated temperature dependence of volume in LTA and CHA. Vibrational density of states of LTA. 6. Acknowledgments The authors declare no competing financial interest.



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 21

References 1. Kärger, J.; Ruthven, D. M.; Theodorou, D. N., Eight-Ring Zeolites. In Diffusion in Nanoporous Materials, Wiley-VCH Verlag GmbH & Co. KGaA: 2012; pp 561-606. 2. Krishna, R.; van Baten, J. M., A Molecular Dynamics Investigation of the Diffusion Characteristics of Cavity-Type Zeolites with 8-Ring Windows. Micropor Mesopor Mat 2011, 137, 83-91. 3. Davis, M. E., Ordered Porous Materials for Emerging Applications. Nature 2002, 417, 813-821. 4. Wang, Z. P.; Yu, J. H.; Xu, R. R., Needs and Trends in Rational Synthesis of Zeolitic Materials (Vol 41, Pg 1729, 2012). Chem Soc Rev 2012, 41, 8210-8210. 5. Van Speybroeck, V.; Hemelsoet, K.; Joos, L.; Waroquier, M.; Bell, R. G.; Catlow, C. R. A., Advances in Theory and Their Application within the Field of Zeolite Chemistry. Chem Soc Rev 2015, 44, 7044-7111. 6. Palomino, M.; Cantin, A.; Corma, A.; Leiva, S.; Rey, F.; Valencia, S., Pure Silica Itq32 Zeolite Allows Separation of Linear Olefins from Paraffins. Chem Commun 2007, 12331235. 7. Krishna, R.; van Baten, J. M., A Molecular Dynamic Investigation of the Diffusion of Methane-Ethane and Methane-Propane Mixtures in Zeolites. Chem Eng Technol 2006, 29, 1429-1437. 8. Hedin, N.; DeMartin, G. J.; Roth, W. J.; Strohmaier, K. G.; Reyes, S. C., Pfg Nmr Self-Diffusion of Small Hydrocarbons in High Silica Ddr, Cha and Lta Structures. Micropor Mesopor Mat 2008, 109, 327-334. 9. Hedin, N.; DeMartin, G. J.; Strohmaier, K. G.; Reyes, S. C., Pfg Nmr Self-Diffusion of Propylene in Itq-29, Caa and Nacaa: Window Size and Cation Effects. Micropor Mesopor Mat 2007, 98, 182-188. 10. Selassie, D.; Davis, D.; Dahlin, J.; Feise, E.; Haman, G.; Sholl, D. S.; Kohen, D., Atomistic Simulations of Co2 and N-2 Diffusion in Silica Zeolites: The Impact of Pore Size and Shape. J Phys Chem C 2008, 112, 16521-16531. 11. Awati, R. V.; Ravikovitch, P. I.; Sholl, D. S., Efficient and Accurate Methods for Characterizing Effects of Framework Flexibility on Molecular Diffusion in Zeolites: Ch4 Diffusion in Eight Member Ring Zeolites. J Phys Chem C 2013, 117, 13462-13473. 12. Boulfelfel, S. E.; Ravikovitch, P. I.; Sholl, D. S., Modeling Diffusion of Linear Hydrocarbons in Silica Zeolite Lta Using Transition Path Sampling. J Phys Chem C 2015, 119, 15643-15653. 13. Combariza, A. F.; Gomez, D. A.; Sastre, G., Simulating the Properties of Small Pore Silica Zeolites Using Interatomic Potentials. Chem Soc Rev 2013, 42, 114-127. 14. Bueno-Perez, R.; Calero, S.; Dubbeldam, D.; Ania, C. O.; Parra, J. B.; Zaderenko, A. P.; Merkling, P. J., Zeolite Force Fields and Experimental Siliceous Frameworks in a Comparative Infrared Study. J Phys Chem C 2012, 116, 25797-25805. 15. Shi, H. B.; Migues, A. N.; Auerbach, S. M., Ab Initio and Classical Simulations of the Temperature Dependence of Zeolite Pore Sizes. Green Chem 2014, 16, 875-884. 16. Sahoo, S. K.; Nair, N. N., A Potential with Low Point Charges for Pure Siliceous Zeolites. J Comput Chem 2015, 36, 1562-1567. 17. Hill, J. R.; Sauer, J., Molecular Mechanics Potential for Silica and Zeolite Catalysts Based on Ab-Initio Calculations .1. Dense and Microporous Silica. J Phys Chem-Us 1994, 98, 1238-1244. 18. Hill, J. R.; Sauer, J., Molecular Mechanics Potential for Silica and Zeolite Catalysts Based on Ab-Initio Calculations .2. Aluminosilicates. J Phys Chem-Us 1995, 99, 9536-9550.


Page 19 of 21

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


19. Sanders, M. J.; Leslie, M.; Catlow, C. R. A., Interatomic Potentials for Sio2. J Chem Soc Chem Comm 1984, 1271-1273. 20. Schroder, K. P.; Sauer, J., Potential Functions for Silica and Zeolite Catalysts Based on Ab Initio Calculations .3. A Shell Model Ion Pair Potential for Silica and Aluminosilicates. J Phys Chem-Us 1996, 100, 11043-11049. 21. Nasluzov, V. A.; Ivanova, E. A.; Shor, A. M.; Vayssilov, G. N.; Birkenheuer, U.; Rosch, N., Elastic Polarizable Environment Cluster Embedding Approach for Covalent Oxides and Zeolites Based on a Density Functional Method. J Phys Chem B 2003, 107, 22282241. 22. Vanbeest, B. W. H.; Kramer, G. J.; Vansanten, R. A., Force-Fields for Silicas and Aluminophosphates Based on Abinitio Calculations. Phys Rev Lett 1990, 64, 1955-1958. 23. Ramsahye, N. A.; Bell, R. G., Cation Mobility and the Sorption of Chloroform in Zeolite Nay: Molecular Dynamics Study. J Phys Chem B 2005, 109, 4738-4747. 24. Demontis, P.; Suffritti, G. B.; Quartieri, S.; Fois, E. S.; Gamba, A., MolecularDynamics Studies on Zeolites .3. Dehydrated Zeolite-A. J Phys Chem-Us 1988, 92, 867-871. 25. Nicholas, J. B.; Hopfinger, A. J.; Trouw, F. R.; Iton, L. E., Molecular Modeling of Zeolite Structure .2. Structure and Dynamics of Silica Sodalite and Silicate Force-Field. J Am Chem Soc 1991, 113, 4792-4800. 26. Krokidas, P. G.; Skouras, E. D.; Nikolakis, V.; Burganos, V. N., Lattice Dynamics Simulation of Thermal Contraction of Faujasites. J Phys Chem C 2010, 114, 22441-22448. 27. Lopes, P. E. M.; Murashov, V.; Tazi, M.; Demchuk, E.; MacKerell, A. D., Development of an Empirical Force Field for Silica. Application to the Quartz-Water Interface. J Phys Chem B 2006, 110, 2782-2792. 28. Zimmerman, P. M.; Head-Gordon, M.; Bell, A. T., Selection and Validation of Charge and Lennard-Jones Parameters for Qm/Mm Simulations of Hydrocarbon Interactions with Zeolites. J Chem Theory Comput 2011, 7, 1695-1703. 29. Garcia-Sanchez, A.; Dubbeldam, D.; Calero, S., Modeling Adsorption and SelfDiffusion of Methane in Lta Zeolites: The Influence of Framework Flexibility. J Phys Chem C 2010, 114, 15068-15074. 30. Krishna, R.; van Baten, J. M., Comment on "Modeling Adsorption and Self-Diffusion of Methane in Lta Zeolites: The Influence of Framework Flexibility". J Phys Chem C 2010, 114, 18017-18021. 31. Bornhauser, P.; Bougeard, D., Intensities of the Vibrational Spectra of Siliceous Zeolites by Molecular Dynamics Calculations. I. Infrared Spectra. J Phys Chem B 2001, 105, 36-41. 32. Corma, A.; Rey, F.; Rius, J.; Sabater, M. J.; Valencia, S., Supramolecular SelfAssembled Molecules as Organic Directing Agent for Synthesis of Zeolites. Nature 2004, 431, 287-290. 33. Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J., Cp2k: Atomistic Simulations of Condensed Matter Systems. Wires Comput Mol Sci 2014, 4, 15-25. 34. VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J., Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput Phys Commun 2005, 167, 103-128. 35. VandeVondele, J.; Hutter, J., An Efficient Orbital Transformation Method for Electronic Structure Calculations. J Chem Phys 2003, 118, 4365-4369. 36. Krack, M.; Parrinello, M., All-Electron Ab-Initio Molecular Dynamics. Phys Chem Chem Phys 2000, 2, 2105-2112. 37. Lippert, G.; Hutter, J.; Parrinello, M., The Gaussian and Augmented-Plane-Wave Density Functional Method for Ab Initio Molecular Dynamics Simulations. Theor Chem Acc 1999, 103, 124-140. 19 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 21

38. Goedecker, S.; Teter, M.; Hutter, J., Separable Dual-Space Gaussian Pseudopotentials. Phys Rev B 1996, 54, 1703-1710. 39. Hartwigsen, C.; Goedecker, S.; Hutter, J., Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn. Phys Rev B 1998, 58, 3641-3662. 40. Becke, A. D., Density-Functional Exchange-Energy Approximation with Correct Asymptotic-Behavior. Phys Rev A 1988, 38, 3098-3100. 41. Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H., Results Obtained with the CorrelationEnergy Density Functionals of Becke and Lee, Yang and Parr. Chem Phys Lett 1989, 157, 200-206. 42. Lee, C. T.; Yang, W. T.; Parr, R. G., Development of the Colle-Salvetti CorrelationEnergy Formula into a Functional of the Electron-Density. Phys Rev B 1988, 37, 785-789. 43. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalized Gradient Approximation Made Simple. Phys Rev Lett 1996, 77, 3865-3868. 44. Frenkel, D.; Smit, B., Understanding Molecular Simulation : From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, 2002, p xxii, 638 p. 45. Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular Dynamics. Journal of Computational Physics 1995, 117, 1-19. 46. Baerlocher, C.; McCusker, L. B. Database of Zeolite Structures. http://www.izastructure.org/databases/. 47. Villaescusa, L. A.; Bull, I.; Wheatley, P. S.; Lightfoot, P.; Morris, R. E., The Location of Fluoride and Organic Guests in 'as-Made' Pure Silica Zeolites Fer and Cha. J Mater Chem 2003, 13, 1978-1982. 48. Martinez-Inesta, M. M.; Lobo, R. F., Investigation of the Negative Thermal Expansion Mechanism of Zeolite Chabazite Using the Pair Distribution Function Method. J Phys Chem B 2005, 109, 9389-9396. 49. Woodcock, D. A.; Lightfoot, P.; Villaescusa, L. A.; Diaz-Cabanas, M. J.; Camblor, M. A.; Engberg, D., Negative Thermal Expansion in the Siliceous Zeolites Chabazite and Itq4: A Neutron Powder Diffraction Study. Chem Mater 1999, 11, 2508-2514. 50. Carey, T.; Corma, A.; Rey, F.; Tang, C. C.; Hriljac, J. A.; Anderson, P. A., The Effect of Extra Framework Species on the Intrinsic Negative Thermal Expansion Property of Zeolites with the Lta Topology. Chem Commun 2012, 48, 5829-5831. 51. Carey, T.; Tang, C. C.; Hriljac, J. A.; Anderson, P. A., Chemical Control of Thermal Expansion in Cation-Exchanged Zeolite A. Chem Mater 2014, 26, 1561-1566. 52. Wragg, D. S.; Morris, R. E.; Burton, A. W., Pure Silica Zeolite-Type Frameworks: A Structural Analysis. Chem Mater 2008, 20, 1561-1570. 53. Verploegh, R. J.; Nair, S.; Sholl, D. S., Temperature and Loading-Dependent Diffusion of Light Hydrocarbons in Zif-8 as Predicted through Fully Flexible Molecular Simulations. J Am Chem Soc 2015, 137, 15760-15771. 54. Awati, R. V.; Ravikovitch, P. I.; Sholl, D. S., Efficient Calculation of Gas Diffusivity in Single-Component and Binary Mixtures of Spherical Adsorbates in Flexible 8mr Zeolites. J Phys Chem C 2015, 119, 16596-16605. 55. Camp, J. S.; Sholl, D. S., Transition State Theory Methods to Measure Diffusion in Flexible Nanoporous Materials: Application to a Porous Organic Cage Crystal. The Journal of Physical Chemistry C 2016, 120, 1110-1120.


Page 21 of 21

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


TOC Graphic


Improved HillâSauer Force Field for Accurate ... - ACS Publications

Recommend Documents

Improved HillâSauer Force Field for Accurate ... - ACS Publications