Molecular Dynamics Study of Diffusion and Surface Permeation of

Subscriber access provided by - Access paid by the | UCSB Libraries

C: Surfaces, Interfaces, Porous Materials, and Catalysis

A Molecular Dynamics Study of Diffusion and Surface Permeation of Benzene in Silicalite German Sastre, Jörg Kärger, and Douglas M. Ruthven J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b00520 • Publication Date (Web): 09 Mar 2018 Downloaded from http://pubs.acs.org on March 9, 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.

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

A Molecular Dynamics Study of Diffusion and Surface Permeation of Benzene in Silicalite

German Sastrea, Jörg Kärger b*, Douglas M. Ruthven c*

a

Instituto de Tecnologia Quimica U.P.V.-C.S.I.C. Universidad Politecnica de Valencia. Avenida Los Naranjos s/n, 46022 Valencia (Spain). b

Faculty of Physics and Earth Sciences, University of Leipzig, Linnéstraße 5, 04103 Leipzig, Germany. E-mail: [email protected] c

Department of Chemical and Biological Engineering, University of Maine, Orono, ME, USA. E-mail: [email protected]

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

Abstract Molecular dynamic simulations have been carried out to determine uptake and release rates for benzene in an idealized crystal of silicalite (the pure silica form of ZSM-5) with two external surfaces perpendicular to the straight channel [010]. A realistic model has been developed in order to simulate a system in which the kinetics are controlled by the combined effects of surface resistance (pore blocking at the external surface) and intra-crystalline diffusional resistance. The system has been treated using periodic boundary conditions and contains a finite reservoir in which (for uptake calculations) the benzene molecules are initially located. For the calculation of release, the benzene molecules are all located initially within the crystal, and the periodicity along the length of the reservoir is removed so that molecules are released at zero pressure. The effect of a surface barrier has been investigated by considering three systems with different degrees of channel entrance blocking (0%, 50% and 87.5%). The resulting calculations of uptake and release make it possible to estimate the relative importance of surface resistance (pore blocking) and intra-crystalline diffusion in determining the sorption rate. It is shown that the classical model based on the formal solution of the one-dimensional diffusion equation, taking account of the finite rate of permeation through the crystal surface (i.e. surface resistance) via the boundary condition, provides a good representation of the kinetic behaviour. For comparison self-diffusion and tracer exchange are also simulated for the same system under comparable conditions.

2 Environment ACS Paragon Plus

Page 2 of 25

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


1. Introduction. Nanoporous materials have found widespread application in catalysis, separation processes, drug delivery, energy storage and gas sensing devices. Control of the size and shape of the nanoporous particles and their transport properties is critical to maximising the efficiency of such processes. In these systems transport of molecules into and out of the nanopores is often the rate-limiting step so, from both fundamental and practical perspectives, it is clearly important to understand the dynamic behaviour at the molecular scale. Recent studies1 have revealed that sorption rates often show substantial variation between apparently identical crystals within the same sample, thus highlighting the need for detailed characterization. It is often assumed that the sorption rate is controlled entirely by diffusion into or out of the particle, with equilibrium between the surface layer and the surrounding fluid. For many systems this is a reasonable approximation. However, following the message of early pulsed field gradient nuclear magnetic resonance (PFG NMR) studies2, recent more detailed studies by advanced optical techniques such as interference microscopy (IFM) and infrared microscopy (IRM)3

4 5

have revealed that, in many zeolite crystals, there is significant mass

transfer resistance at the external surface so that the sorption rate is actually controlled by the combined effects of internal diffusion and surface resistance. This has also been confirmed by macroscopic techniques such as ZLC (zero length chromatography)6 7. Variations in surface and internal resistance can be explained by pore blocking and structural defects. A new picture has emerged for real zeolite crystals for which, in some cases, transport rates over long distances (many unit cells) reflect the influence of surface and internal barriers rather than the pore structure of the idealized framework. As a result, the apparent intra-crystalline diffusivities often show a strong dependence on the length scale of the measurement, as first observed by PFG NMR8, modelled with Monte Carlo techniques9, and confirmed by combined transmission electron microscopy, PFG NMR and QENS (quasielastic neutron scattering) 10. Confocal fluorescence microscopy has revealed a large range of different sub-structures within what appear to be well formed MFI-type zeolite crystals11

12

but the impact of such structural defects on intra-crystalline diffusion has not yet been determined.



Page 4 of 25

In the present study, we model the behaviour of three ideal crystals of silicalite bounded by two external surfaces with different degrees of pore blocking, in order to study how the surface resistance affects the overall sorption kinetics.

2. Methodology and Models. 2.1. Structure of Silicalite. We selected the benzene–silicalite system because of the large amount of accumulated experimental data13

14 15 16 17 18 19

and many computer simulations20

21 22 23

that are available

for comparison with our results. Three different structures of silicalite have been considered according to published material24

25

26,

with respective space groups: Pnma (#62,

orthorhombic, also from ref. 27), P212121 (#19, orthorhombic) and P21/c (#14, monoclinic), of which the second (so called PARA-silicalite)17 was selected, based on a previous study28 in which the reproducibility of pore sizes was modelled using different force fields. This structure was solved taking into account the effect of p-xylene loading on the structure, so the structural model accounts for the deformation of the channels resulting from the presence of the guest molecules. 2.2. Calculation of 10-ring Diameters. Silicalite contains two types of channels (straight and sinusoidal). In the case of PARA-silicalite there are four different 10-rings in the sinusoidal channel (along [100]) and three different 10-rings along the straight channel (along [010]). Each type of ring is analysed independently to ensure a match between the calculated and experimental dimensions. A new version (1.2.7.6) of the zeoTsites29 software was employed for the analysis and classifications of all 10-rings in the ZSM-5 model28. The 10-ring diameters are calculated dynamically, for each saved conformation of the molecular dynamics runs (vide infra), so that the effect of thermal vibrations on the calculated ring diameters is properly accounted for. This should provide a more reliable comparison with the XRD reported ring diameters (sinusoidal 5.1×5.5 Å, and straight 5.3×5.6 Å) whose values are also thermally averaged. 2.3. Molecular Dynamics. A 2×2×2 unit cell (a = 2×20.121 Å, b = 2×19.82 Å, c = 2×13.438 Å) of silicalite (containing 8×96 SiO2 atoms) was used as the model. The cell length along [010] was enlarged to 79.46 Å and the resulting two surfaces perpendicular to [010] were terminated by


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


silanol groups (≡SiOH). This model (Figure 1) contains 8 straight channels, each with two external surfaces, thus offering possible degrees of surface blocking between 0 and 16. The extended length (79.46 Å) along [010] corresponds to 43.8 Å of zeolite crystal (with silanol groups at both ends) and 35.66 Å of reservoir. The periodic nature of this crystal means that the two ends of the reservoir are joined and the periodic model is a repetition of [crystalreservoir] units.

Figure 1. Model of a silicalite crystal employed for the calculations. Blocking atoms are highlighted. The half-length within the crystal (between the reservoirs) is ℓ = 21.9 Å.

In the simulations of uptake the reservoir is initially filled with a given number of thermalized benzene molecules while the zeolite crystal contains no benzene. The simulation starts with a random configuration of benzene molecules, all of them located in the reservoir, and from this configuration the system is allowed to evolve by itself (at constant volume, 298 K, and with a fixed number of particles), without any restraint or external force. As a consequence of the lower free energy of the benzene-zeolite (intracrystalline) with respect to the benzene-benzene (reservoir) systems, benzene molecules tend to adsorb and diffuse inside the zeolite micropores, until the external pressure becomes low enough that it approaches the equilibrium state. When simulating release the starting configuration is the final equilibrated configuration from the previous (adsorption) run, but with all molecules in the reservoir removed. Another important aspect of release simulation is the removal of periodicity along



Page 6 of 25

the [010] direction, leading to a system where the released molecules escape into an infinite reservoir maintained at zero pressure. DL_POLY 2.2030 was used for all the molecular dynamic simulations reported in this work including full flexibility and periodic boundary conditions for all the atoms of the system, except the blocking atoms which remain in the same position at the centre of the channel (just inside the external surface). The temperature chosen is 298 K within the NVT ensemble. We employed the Verlet-leapfrog integration algorithm and the Evans thermostat, with a timestep of 1 fs. Each run comprised an equilibration stage of 5×104 steps followed by the necessary production stage so as to ensure sufficiently good statistics for the analysis, which amounts to simulations times between 100-1000 ns depending on the run. Since all benzene molecules are initially located in the reservoir, the equilibration of their velocities is very fast; less than 50 ps in all cases. The length of the equilibration does not need to increase accordingly to the total time of the simulation and hence 50 ps is sufficient regardless the total length of the simulation. The cut-off for the non-bonding forces was set to 9 Å, and the Ewald summation was employed for the Coulombic with a precision parameter set to 10 -3. The configurations were saved every 1000 time-steps (1 ps). 2.4. Force Field . Five full-atom force fields were tested previously to benchmark the accuracy of dynamic pore size in ZSM-528. The results obtained with the force fields by van Best, Kramer and van Santen31 (BKS); Nicholas et. al.32 (Nic); Sastre33

34

(Sas); Vessal, Leslie, Catlow35

(VLC); and Sanders, Leslie, Catlow36 (SLC), suggest that Sas force field, (with charges of 2.1 and -1.05 for Si and O bulk, and, 0.425 and -0.95 for the terminal H and O in the silanol groups), is a good compromise between accuracy and computational efficiency. This was therefore chosen. Although SLC performs better, it takes ca. 5 times longer than the others28. The intramolecular all-atom fully flexible force field by Oie et al.37 has been considered for the benzene molecules, with charges of -0.11 and 0.11 for the C and H atoms respectively. For the benzene-benzene intermolecular interactions the all-atom OPLS 38 was used, which has been specifically parameterised using data from liquid benzene. For the zeolite-benzene interactions we used the force field developed by Snurr et al.39, which was recently found the best choice in a study of diffusion of benzene in ZSM-5 22. The blocking atoms have been fixed at the centre of those 10-ring straight channels where they are located, at a distance from the external surface slightly larger than 4 Å (Figure


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


1). A Lennard-Jones force field corresponding to a large atom has been found convenient so that the benzene molecules cannot penetrate the channel space between the atom and the pore wall, and to this effect the parameters for the Rb atom from the Universal force field40 have been selected. The interactions of these blocking atoms (of neutral charge) with the C and H atoms of benzene were obtained using the Lorentz-Berthelot mixing rules. 2.5. Calculation of Uptake and Release. Three regions have been considered in our model (Figure 1); the reservoir, the external surface of the silicalite crystal and the bulk crystal. Only the molecules that have penetrated deeper than 4 Å from the outer H atoms (belonging to silanol groups) at the external surface, are considered to be inside the crystal. The rest of the molecules i.e. those which are located within the surface region, plus the molecules in the reservoir (with which they are in rapid dynamic equilibrium) are considered to be outside the crystal. At any time during the molecular-dynamics run it is possible to calculate the number of molecules inside/outside, giving the uptake or release curves for each simulation. Intra-crystalline self-diffusion coefficients are calculated from the mean square displacement (MSD) using the equations from previous work28, and considering a time span corresponding to the equilibrium region of the uptake curve. Only the benzene molecules which remain inside the crystal over the selected time period are considered. 2.6. Modelling and Quantification of Transient Sorption and Exchange Rates. To quantify the sorption and exchange rates generated by the simulations the transient uptake or exchange curves were matched to the relevant solutions of the one-dimensional diffusion equation for the appropriate initial and boundary conditions as summarized by Crank41. The relevant solutions are given in Crank’s Figures 4.6 (for a diffusion controlled system) and 4.7 (for a system controlled by the combined effect of internal diffusion and surface resistance). These solutions assume a constant diffusivity (independent of concentration) but, except at low concentrations, one may expect the diffusivity in an adsorbed phase to be concentration dependent. However, the shape of the initial portion of the transient uptake curve (up to about 50% fractional approach to equilibrium) is relatively insensitive to concentration dependence of the diffusivity. Thus the constant diffusivity model can be used as a reasonable approximation even when the diffusivity is concentration dependent. However, the resulting integral diffusivity will be a weighted average of the diffusivity over the relevant concentration range (see ref.42 Section 11.7). Alternatively, where



Page 8 of 25

the uptake or release curve is available over the full range, the time constant (and hence the effective diffusivity) may be estimated by the method of moments. The values obtained by both approaches are reasonably consistent (see Supporting Information, section S6). Several different diffusivities (transport diffusivity, tracer diffusivity, self-diffusivity) are derived from the simulations depending on the nature of the system. A short summary of the definitions and the relationship between these parameters is given in the Supporting Information (section S1). The validity of a diffusion model depends on statistical averaging which allows molecular concentrations and fluxes to be treated as continuous variables. For a large system containing many molecules this is essentially exact. However, the systems considered in the present study are relatively small (eight unit cells containing only between 40 and 64 molecules) so exact statistical averaging and therefore exact conformity with a diffusion model cannot be expected. Nevertheless the diffusion models are found to provide a good qualitative (and approximately quantitative) representation of the observed behaviour. The parameters derived by matching the uptake/exchange curves to Crank’s equations are summarized in Table 1. Table 1. Parameters derived from the fit of Crank’s equations. Figure 2/S2

Loading (molec./u.c.) adsorption 7.5 Mode

3/S4,S9 adsorption 6 5(a) 5(b)

desorption exchange exchange

Λ 0.19

5.9

0.78

8 5.8 5.8

0 0.73 0

% block 0 50 87.5 0 50 87.5 0 0 0

D/ℓ2 (s-1)

L=kℓ/D

k (cm·s-1)

3.75×106 3.75×106 3.75×106 1.44×106 1.44×106 1.44×106 1.25×106 6.25×106 6.25×106

∞ 4 1 ∞ 0.5 0.125 ∞ ∞ ∞

3.3 0.82 0.16 0.04 -

3. Results and Discussion. 3.1. Uptake of 310 Benzene Molecules in Silicalite. The choice of 310 benzene molecules was made in order to simulate a density similar to liquid benzene in the external fluid at equilibrium. Three different degrees of blocking have been considered with 0, 8 and 14 blocking atoms in a total of 16 openings to the straight channel along [010], which corresponds to fractional blocking of 0 %, 50 % and 87.5 % of the pore entrances. Starting from 310 benzene molecules in the reservoir while the crystal is empty, the uptake curves, generated in a MD simulation run of 120 ns, are shown in Figure 2.


Page 9 of 25 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 three systems reach similar equilibrium uptakes, with N∞ ≈ 60 benzene molecules inside the crystal, containing 2×2×2 unit cells (7.5 benzene/u.c.), although with different kinetics. This is close to the experimental saturation capacity for benzene reported by Olson et al. (7.6 molecules/unit cell)27. Others5,8,13 have reported slightly larger saturation values approaching the theoretical limit of 8 molecules per unit cell. With 60 benzene molecules inside the crystal, the resulting number of molecules in the reservoir will be 250. The reservoir volume (40.242×35.66×26.876 Å3) gives a density of benzene of 841 kg/m3, close to the experimental value of 876 kg/m3 for liquid benzene at 298 K. This indicates that, with the silicalite crystal element shown in Figure 1, 310 molecules is a good choice for the total number of benzene molecules in order to simulate a realistic system, since it yields an equilibrium fluid density in the reservoir close to the actual density of liquid benzene at 298 K. For the unblocked system we assume pure diffusion control. Taking N∞ = 60 (7.5 molecules per unit cell) the fraction of molecules eventually adsorbed from the reservoir is 60/310 = 0.19 = Λ. We therefore estimate the diffusivity from Crank’s Figure 4.6 with Λ = 0.19 (interpolated) and ℓ = 21.9 Å (see caption of Figure 1). Hence D ≈ 1.8×10-7 cm2s-1. The theoretical curve calculated using these parameters provides a good representation of the MD data (except at very short times) as may be seen from Figure 2.

Figure 2. Number of benzene molecules adsorbed (Nt/N∞) in the silicalite crystal as a function of time, with different degrees of pore blocking; unblocked (0% blocked), 50% blocked, 87.5% blocked. The total number of benzene molecules is 310. The lines show the theoretical solutions of Crank’s



equations (Figures 4.6 and 4.7 in ref.41), calculated from the one-dimensional diffusion equation with ℓ = 21.9 Å, D = 1.8×10-7 cm2s-1 and Λ = 0.19 (N∞/Ntot = 60/310) for L = ∞, L = 4, L = 1.

To interpret the data for 50% and 87.5% blocking we use Crank’s Figure 4.7, which summarises the solution for a system in which the uptake rate is controlled by the combined effects of intra-crystalline diffusion and surface resistance (ignoring the effect of the small time variation in the boundary condition). Since the intra-crystalline diffusivity should not be affected by pore blocking we assume D = 1.8×10-7 cm2s-1 as for the unblocked crystal. An iterative procedure was then used to estimate the values of the parameter L = kℓ/D which best represent the MD uptake curves. Thus we find: L ≈ 4, k ≈ 3.3 cm·s-1 (50% blocking), and L ≈ 1, k ≈ 0.82 cm·s-1 (87.5% blocking). We note that these values appear to be consistent since k scales, as expected, with the fraction of unblocked channels (i.e. with 1 - Fractional Blocking, see Supporting Information S7). The theoretical uptake curves re-calculated with these parameters provide reasonable representations of the MD simulation data, as may also be seen from Figure 2. There is some deviation at short times, especially for 50% blocking, but the MD data for this run show an almost immediate jump to 10 molecules adsorbed at very short times, which appears to be a statistical anomaly. 3.2. Uptake of 60 Benzene Molecules in Silicalite. A lower loading of 60 benzene molecules in the 2×2×2 unit cell of silicalite (Figure 1) was considered in order to assess the effect of concentration on the sorption kinetics. Initially all 60 benzene molecules are in the reservoir while the crystal is empty. The uptake results are shown in Figure 3. For this system N∞ = 47 and Λ = 47/60 = 0.78 and the equilibrium loading corresponds to 5.9 molecules per unit cell, a fractional loading of 74% (based on the theoretical saturation limit of 8 mol./u.c.). In principle it should be possible to estimate the equilibrium benzene pressure from the long time asymptote of the uptake curve. However, for the present system this is not practical owing to the small size of the model system and the strong adsorption of benzene. This is discussed in the Supporting Information (section S11). The kinetic data for the unblocked system are analysed using Crank’s Figure 4.6 with Λ = 0.78 and ℓ = 21.9 Å (Supporting Information, S4). This yields D ≈ 6.9×10-8 cm2s-1, roughly half of the value obtained (1.8×10-7 cm2s-1) at the higher loading (7.5 molec./u.c. or 94% fractional loading), suggesting that the diffusivity increases with loading, as expected. The theoretical curve re-calculated using these parameters provides a reasonable representation of the MD data (Figure 3).


Page 10 of 25

Page 11 of 25 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 high value of Λ makes it difficult to analyse the uptake curves for the blocked systems in detail. For pure diffusion control (zero blocking), since the value of Λ is known, it is possible to estimate (from Crank’s Figure 4.6 or Eq. 4.37 in ref. 41) the value of D/ℓ2 (and hence D) that provides a good fit of the data. When Λ is small, using Crank’s Figure 4.7 (Eq. 4.53), we can estimate the value of L (kℓ/D), and hence k, since we know D from the zero blocking data. The difficulty arises when Λ is large since we do not have the solution for combined diffusion and surface resistance except for Λ = 0. To solve this problem in an approximate way one possibility is to use the correction factor from moments analysis which suggests t/to = (1 – Λ), where t is the time for a given fractional uptake and to is the time for the same fractional uptake when Λ = 0. However, although this is reasonable when Λ is small, it breaks down when Λ is large, as it is for some of our systems. Instead, we therefore used Crank’s Figure 4.6 to estimate the ratio t/to as a function of the fractional approach to equilibrium. Assuming that the same factor applies when the rate is controlled by the combined effects of diffusion and surface resistance, we can use these factors to estimate the dimensionless time, and hence L values for the blocked samples, from Crank’s Figure 4.7. A detailed calculation is given in Supporting Information (section S9). For the 50% blocked sample, if we take Λ = 0.78, L = kℓ/D = 0.5 with D = 6.9×10-8 cm2s-1 and ℓ = 21.9 Ǻ we get k = 0.16 cm∙s-1 which is much smaller than the value (k = 3.3 cm∙s-1) obtained for the higher loading case (Figure 2).

Figure 3. Number of benzene molecules adsorbed in the silicalite crystal (Nt/N∞) as a function of square root of time, with different degrees of pore blocking; unblocked (0% blocked), 50% blocked, 87.5% blocked. The total number of benzene molecules is 60 and the equilibrium loading is 5.9



Page 12 of 25

molec./u.c.. The broken line for 0% block shows the theoretical solution of Crank’s equation (Figure 4.6 in ref. 41) calculated from the one-dimensional diffusion equation with ℓ = 21.9 Å, D = 6.9×10-8 cm2s-1 and Λ = 0.78 (N∞/Ntot = 47/60). The dotted line for 50% block shows the theoretical solution of Crank’s equation (Figure 4.7 in ref.41), with the above ℓ and D values and L = 0.5.

Since 1/k represents the increase in diffusion time required to circumnavigate the blockage at the pore entrance k should be proportional to the differential transport diffusivity close to the external surface. In an adsorption step the loading at the external surface will correspond, at all times, to the equilibrium loading. If D increases with loading, as the comparison between Figures 2 and 3 suggests, we should expect k to be greater for the higher loading case. The loading dependence of the differential diffusivity will be much stronger than that of the integral diffusivity. The ratio 3.3/0.158 = 21 (going from 94% to 74% fractional loading) seems rather large but still within the physically reasonable range. If L= 0.5 for 50% blocking we would expect L ≈ 0.125 for 87.5% blocking. At this small L value the uptake rate will be controlled almost entirely by surface resistance, so the uptake curve should approach the limiting rate expression: Nt/N∞ = 1 – exp[-kt/(1 – Λ)ℓ] = 1 – exp[-3.28×106t(s)]

(1)

and at short times: Nt/N∞ ≈ 3.28×106t(s)

(2)

This expression correctly predicts the initial uptake rate but not the shape of the uptake curve. Due to the very low uptake the simulation data for this case are not statistically reliable. 3.3. Self-Diffusivities (D ). The self-diffusivities may be conveniently estimated from the slopes of plots of MSD vs time: d(MSD)/dt = 6 D ; d(MSD(x))/dt = 2 Dx ; d(MSD(y))/dt = 2 Dy. Such plots are shown in the Supporting Information (S3 and S5). The resulting self-diffusivities for both loading levels (7.5 and 5.9 molecules per unit cell) are summarised in Table 2. The MSD is always substantially smaller than ℓ2 so the surface blocking should have no effect on the selfdiffusivity. In all but one case the slopes of the plots of MSD(x) and MSD(y) vs time are essentially the same showing that Dx ≈ Dy as has been found experimentally for this system19. In general the values follow the expected trend: Dx ≈ Dy> D >> Dz.


Page 13 of 25 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.4. Tracer Exchange. During the equilibrium period (Figure 4), provided that there is sufficient exchange of molecules between the intra-crystalline region and the surrounding fluid, one may treat the exchange process as a “tracer exchange” and derive the tracer diffusivity (which should correspond to the self-diffusivity) by matching the exchange curve from the MD simulation to the appropriate solution of the one-dimensional diffusion equation. In contrast to the transient uptake curves, the total loading is fixed so the diffusivity will be constant and should correspond to the self-diffusivity at the same loading. The required conditions are fulfilled for the unblocked crystal with Ntot = 60, N∞ = 44.

Figure 4. Uptake of benzene in silicalite after equilibrium is reached as shown by the approximately stationary uptake. The total number of benzene molecules is 60.

Two different simulations (Figure 5) were performed at the same equilibrium loading (Ntot = 60, N∞ = 44 and 5.5 molecules/unit cell). In simulation (a), which mimics a normal tracer exchange experiment, a molecule initially inside the crystal is considered to be “desorbed” as long as it remains outside the crystal but if it re-enters the crystal it is redesignated as an “adsorbed” molecule (and vice versa for the molecules initially outside the crystal). In simulation (b) a molecule that leaves the crystal is considered as desorbed even if it re-enters the crystal some time later. In this way the surface concentration of marked molecules (those which were inside the crystal at time zero) is held artificially at zero while the total number of molecules remains constant. This mimics a very rapid irreversible reaction at the surface or, equivalently, desorption into an infinitely large reservoir. Since these two



situations [(a) and (b)] differ only in the boundary conditions at the external surface they should be described by the same diffusivity (with properly chosen boundary conditions). The (a) simulations proved straightforward. However the application of the (b) boundary condition proved more challenging. A molecule arriving at the external surface has a high probability of returning very quickly to the site from which it came. Consequently, if one counts a molecule as desorbed (or adsorbed) from the instant that it first arrives at the surface one obtains an unrealistically high rate of exchange between the phases. This was investigated by varying the time (or equivalently the number of consecutive configurations, x) for which a molecule must remain in the new phase before it is counted as having changed phase. The results are summarised in the Supporting Information (section S10). In time units corresponding to successive configurations it was found that, for x < 5, the sorption rate decreases with increasing x, but for 5 < x > Dx. 43 The thermodynamically corrected (differential) transport diffusivity is generally similar to the self-diffusivity (same order of magnitude). The self-diffusivities derived from the simulations at loadings of 5.5 and 7.5 molecules per unit cell are approximately 1.7×10-8 cm2s-1 and, as expected, this is similar to the tracer exchange diffusivity (3×10-8 cm2s-1 ) - see

Figure 5. The experimental self-

diffusivities derived from tracer exchange and TZLC (tracer zero length chromatography) measurements are smaller by a factor of about two orders of magnitude (≈1-3×10-10cm2s-1). The simulations are based on a perfect crystal whereas real crystals generally contain structural defects, leading to dislocation of the channels with consequent reduction in the diffusivity. The difference between the experimental and simulation values can be explained



on this basis although small errors in the force-field and/or the assumed molecular dimensions may also contribute. The experimental transport diffusivities show a strong and complex concentration dependence and this makes detailed comparisons more difficult. In the low loading limit (the Henry’s law region) the self-diffusivity and transport diffusivity should coincide and the experimental data for benzene-silicalite conform to this expectation. At 298 K the limiting transport diffusivity at low loadings (Do) is similar to the self-diffusivity (≈1-3×10-10cm2s-1) as measured by tracer exchange. Beyond the Henry’s law region the transport diffusivity increases rapidly with loading, mainly as a result of the thermodynamic factor (dlnp/dlnc) which corrects for the non-linearity of the relationship between concentration and thermodynamic activity (as defined by the equilibrium isotherm). However, as a result of molecular re-arrangement13

44 45 46,

the equilibrium isotherm for benzene-silicalite (at 303 K)

shows an inflection at about half saturation (θ = 0.5), leading to a maximum value (≈ 12) for the thermodynamic factor and a corresponding maximum for the transport diffusivity at this loading47

48.

At higher loadings the transport diffusivity again increases as the effect of the

thermodynamic factor becomes dominant. The adsorption simulations (Figures 2 and 3) suggest that, in this region, the integral transport diffusivities also increase with loading: D = 6.9×10-8 cm2s-1 at θ = 0.68 and D = 1.8×10-7 cm2s-1 at θ = 0.94. This is qualitatively consistent with our expectation since the integral diffusivity is directly related to the differential diffusivity at the equilibrium loading. At high loading (θ > 0.9) the integral desorption diffusivity from the MD data (Figure 6) is about one third of the integral adsorption diffusivity (Figure 2) and about double the selfdiffusivity. This is qualitatively as expected for a system in which the transport diffusivity increases strongly with loading, and the ratio Ddes/Dy is approximately consistent with what we would expect (see Fig. 6.17 in ref. 42). It is noteworthy that the adsorption simulations (Figures 2 and 3) show significantly slower uptake rates for the 50% blocked crystal compared with the unblocked crystal whereas the desorption data (Figure 6) show very little difference (although, as expected, the 87.5% blocked crystal is much slower). This difference in behaviour is understandable. The relative importance of surface resistance and intra-crystalline diffusion is determined by the parameter L= kℓ/D, which approaches infinity when there is no surface resistance. The integral diffusivity for desorption is much smaller than that for adsorption with the result that for


Page 18 of 25

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


adsorption L ≈ 4 whereas for desorption L ≈ 14 which is large enough to ensure negligible surface resistance.

4. Conclusions. It has been shown that a small system of 8 unit cells of silicalite, including the crystal surface and a reservoir, can provide a useful model for molecular dynamic (MD) simulations of the kinetic behaviour of adsorbed benzene molecules, allowing the calculation of selfdiffusion and adsorption / desorption kinetics. Despite the numerical difference from the experimental values, self-diffusivities derived from the MD simulations show the same pattern as the experimental data: Dx ≈ Dy> D >> Dz. Integral transport diffusivities derived from adsorption and desorption simulations are much larger than the experimental values, but again show similar trends: Dads > Ddes >> D. Different degrees of pore blocking at the surface have been simulated in order to mimic the observed behaviour of real crystals. Despite the small size of the simulation system and the discrete nature of the surface barrier due to blocking of the pore entrances, the continuum model (internal diffusion plus surface resistance) provides a reasonable representation of the dynamic behaviour. The surface resistance varies with loading in a manner consistent with the concentration dependence of the diffusivity. However, with such a small system, it is only possible to simulate the main features of the kinetic behaviour. To study the effect of structural defects and phenomena such as guest-induced transitions in the host lattice would require a much larger system. Such effects, which have been observed by X-ray diffraction and high resolution NMR spectroscopy as well as from hierarchical atomistic/lattice simulations25

49 50 51,

are beyond the present scope of MD simulations but

present an interesting challenge for the future.

Acknowledgements. G.S. thanks ASIC-UPV computational centre and Paco Rosich for computational assistance, and the Spanish government for the provision of Severo Ochoa (SEV-2016-0683), MAT201571842-P and CTQ2015-70126-R projects.



Supporting Information Available: Definitions of diffusivities, data on mean square displacements, fitting of equations 4.37 and 4.53 from Crank, method of moments to recalculate diffusion coefficients, quantification of surface resistance based on model considerations, replots of Figures 2, 3, 5, 6 in units of time, analysis for a model combining diffusion and surface resistance for Figure 3 (50% blocking), comments on tracer exchange, and adsorption equilibrium for benzene in silicalite.


Page 20 of 25

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


Notation. D

(Transport or Fickian) diffusivity (cm2s-1)

D

Integral diffusivity

D

Self- or tracer diffusivity

k

Surface permeance (cm·s-1), defined by Flux = k·Δc (‘c’ is concentration)

ℓ

Half-length of crystal defined in Figure 1 (ℓ = 21.9 Ǻ)

L

Parameter kℓ/D - ratio of time constants for intra-crystalline diffusion and surface permeation. L>>1 corresponds to diffusion control; L

Molecular Dynamics Study of Diffusion and Surface Permeation of

Recommend Documents