(278a) Computing Diffusivities in Spatially Periodic Media From the Rate Constants of Elementary Jumps Between Sorption Sites: Master Equation Solution by Recursive Reduction of Dimensionality

A quantitative knowledge of the diffusion rates of atoms or molecules within crystalline or amorphous nanoporous or nonporous solids, or on the surfaces of such solids, is essential for the design of sorption separations, heterogeneous catalytic conversions, as well as of doping, morphology control and surface modification processes aimed at improving the properties of materials.  Today, molecular simulation techniques offer a promising avenue for predicting diffusion phenomena based on atomic-level structure and interactions [1].  The most straightforward simulation method for predicting diffusion via molecular simulation is molecular dynamics, in which the trajectories of diffusing atoms or molecules are tracked via numerical solution of the equations of motion for all particles in the system and self- or transport diffusivities are extracted from mean square displacements via appropriate Einstein relations. 

When diffusion is too slow to be tracked by “brute force” molecular dynamics simulation (diffusivity lower than 10^-10 m²/s), methods based on the analysis of infrequent events become valuable.  In many cases diffusive progress is slow because the diffusant spends most of its time trapped within relatively confined sorption sites or “states” within the solid medium and only infrequently jumps from site to site.  Rate constants, i.e. transition probabilities per unit time, can be computed for the elementary jumps, given the potential energy function of the system, via dynamically corrected transition state theory.  Once rate constants have been obtained, the problem of diffusion is reduced to solving the “master equation”,  i.e. a system of first order differential equations in the probabilities of occupancy of the states as functions of time, with prescribed initial and boundary conditions.  In simulation work, periodic boundary conditions are typically used to enable the estimation of macroscopic properties from model systems of length scale 1-10 nm, whether the solid medium under study is crystalline or amorphous.  The connectivity among sites and the rate constants of intersite transitions determine the rate constant matrix appearing in the master equation.  In their 1991 paper. June, Bell, and Theodorou first demonstrated the power of dynamically corrected transition state theory for predicting diffusion in zeolites [2].

The most popular approach for tracking temporal evolution in a network of discrete states connected by interstate transitions with known rate constants is Kinetic Monte Carlo (KMC) simulation.  KMC relies on the generation of a large number of stochastic dynamical trajectories consisting of successive jumps between states [2].  When the spectrum of eigenvalues of the rate constant matrix governing the network is broad, however, KMC becomes inefficient, as it exhausts itself tracking fast transitions and hardly samples slower ones.  Direct solution of the master equation obeyed by the time-dependent probabilities of occupancy of the states is preferable.  In extended spatially periodic systems formed by replicating a “unit cell” of states many times, solution of the master equation calls for diagonalizing a rate constant matrix of  large dimension, which can be challenging computationally. 

This work is dedicated to Professor Alexis T. Bell, a valuable collaborator, mentor, and friend, for his outstanding contributions to the chemical engineering profession.  It takes the ideas presented in reference [2] one step further.  We begin by exploring the structure of the rate constant matrix in a system characterized by spatial periodicity.  We show that the diagonalization of this matrix can be reduced recursively to the diagonalization of a set of matrices of much smaller dimensionality, which refer to a single unit cell in the periodic structure.  Furthermore, we develop a general algorithm for this Master Equation Solution by Recursive Reduction of Dimensionality of the rate constant matrix (MESoRReD). 

To assess the correctness and efficiency of MESoRReD, we apply it to the problem of calculating the low-temperature, low occupancy diffusivity tensor of xenon in the zeolite Silicalite-1, using the states, interstate transitions, and transition state theory-based rate constants computed in Reference [2].  We show that MESoRReD yields correct and accurate diffusivities, while affording savings of more than five orders of magnitude in CPU time relative to KMC.  Furthermore, we demonstrate that, in order to calculate the diffusivity, it suffices to obtain the absolutely smallest nonzero eigenvalue of the rate constant matrix.  In addition, we manage to track which of the small submatrices arising in the recursive reduction scheme contributes this eigenvalue.  With this finding, the problem of calculating diffusivity in a periodic array of states (sites) is ultimately reduced to a computationally trivial determination of the absolutely smallest eigenvalue of a small matrix characterizing a single unit cell, which can be written down by inspection of the transitions occurring within and across the faces of the unit cell. 

Through MeSoRReD, a much-needed step is constructed that leads from topology of a network of states and interstate transition rate constants to macroscopic properties characterizing long-time evolution on the network of states.

References

[1]  J. Kärger, D. M. Ruthven, D. N. Theodorou, Diffusion in Zeolites and Other Microporous Solids, Wiley-VCH, 2012.

[2] R. L. June, A.T. Bell, D.N. Theodorou, J. Phys. Chem., 95, 8866 (1991).