Diffusion of aromatics in silicalite-1: predictions using hierarchical molecular simulations based on transition state theory

Abstract

This doctoral thesis describes the development and application of a new computational methodology for calculating the rate constants governing the motion of sorbate molecules strongly confined within shape-selective nanoporous materials. An umbrella sampling strategy, employing repulsive walls to confine the sorbate within specific regions of the pore space or within specific orientations, is invoked to extract free energy profiles with respect to the sorbate degrees of freedom. Based on these profiles, it is shown how the multidimensional problem of translational diffusion of a sorbate (e.g., benzene) in a flexible nanoporous medium (e.g. in the zeolite silicalite-1) can be reduced first to a 6-dimensional problem, then to a 3-dimensional one, and finally to a 1-dimensional one. A 3-dimensional free energy distribution is accumulated as a function of the sorbate center of mass position and ultimately reduced to a set of 1-dimensional profiles for the sorbate center of mass along the pore axes. From these profiles, the rate constants for jumps executed by the sorbate molecule between sorption sites are calculated using Transition State Theory; from the latter rate constants, the low-occupancy self-diffusivity is obtained. The advantage of this modified methodology of Curvilinear Umbrella Sampling (CUS) in comparison with conventional umbrella sampling, on which it was based, is that it can be applied to curvilinear transition paths in addition to straight ones. In this thesis, a new method was also invented to obtain the low-occupancy self-diffusivity from the rate constants of elementary jumps, in addition to the existing Kinetic Monte Carlo one. This new method solves the master equation for a system evolving on a spatially periodic network of states. The network contains 2ν images of a “unit cell” of n states, arranged along one direction with periodic boundary conditions at the ends. We explore the structure of the symmetrized (2νn)x(2νn) rate constant matrix for this system and derive a recursive scheme for determining its eigenvalues and eigenvectors, and therefore analytically expressing the time-dependent probabilities of all states in the network, based on diagonalizations of nxn matrices formed by consideration of a single unit cell. We call the new method MESoRReD, for Master Equation Solution by Recurive Reduction of Dimensionality. MESoRReD is applied to the problem of low-temperature, low-occupancy diffusion of xenon in the zeolite silicalite-1 using the states, interstate transitions, and transition state theory-based rate constants previously derived by R. L. June et al., J. Phys. Chem., 95, 8866 (1991). MESoRReD yields a diffusion tensor for this system which differs by less than 3% from the values derived previously via Kinetic Monte Carlo (KMC) simulations and confirmed by new KMC simulations conducted in the present thesis. The computational requirements of MESoRReD are compared against those of KMC, numerical solution of the master equation by the Euler method, and direct molecular dynamics simulation. In the problem of diffusion of xenon in silicalite-1, the new method is shown to be faster than these alternative methods by factors of about 3.5x10^5, 1.3x10^6, and 1.7x10^7, respectively. The computational savings and ease of setting up calculations afforded by the MESoRReD in diagonalizing the rate constant matrix make it attractive as a means of predicting medium- and long-time dynamical phenomena in spatially periodic systems from atomic-level information. Next, the specific computational methodologies we have described are applied for the calculation of the self-diffusion coefficients of para-xylene and benzene in silicalite-1 at infinite dilution. In addition, the orientational distributions of phenyl rings and methyl stems of para-xylene and benzene sorbed in the zeolite were studied and the existence of entropic barriers to translational motion was checked. A new reduction method for the states appearing in the free energy profiles is presented and used for the calculation of transition rate constants for elementary jumps. These results agree with Quasi-Elastic neutron scattering measurements. A major conclusion from simulations is that para-xylene diffuses roughly 100 times faster than benzene when sorbed at low occupancy in silicalite. Benzene encounters strong entropic barriers to translational motion at channel intersections, where it can adopt a variety of orientations. The corresponding barriers for para-xylene are much lower, reflecting the inability of its major axis to reorient within channel intersections. Finally, the paradox that the heavier and larger molecule diffuses faster than the lighter and smaller one inside silicalite-1 is explained for first time.

Στην παρούσα διδακτορική διατριβή περιγράφεται η ανάπτυξη και η εφαρμογή μίας νέας μεθοδολογίας για τον υπολογισμό των σταθερών ρυθμού που διέπουν τις κινήσεις των ροφημένων μορίων μέσα σε νανοπορώδη υλικά που λειτουργούν ως μοριακά κόσκινα. Η μεθοδολογία αυτή, βασισμένη στη μεθοδολογία της ‘δειγματοληψίας ομπρέλλας’, εφαρμόζει απωστικούς τοίχους οι οποίοι περιορίζουν το ροφημένο μόριο σε συγκεκριμένες περιοχές των πόρων ενός νανοπορώδους υλικού και βοηθούν στην εξαγωγή των προφίλ της τοπικής ελεύθερης ενέργειας. Με βάση αυτά τα προφίλ, δείχνεται πως το πολυδιάστατο πρόβλημα πρόβλεψης της μεταφορικής διάχυσης του ροφούμενου μορίου μπορεί να αναχθεί πρώτα σ’ένα 6-διάστατο, έπειτα σε ένα 3-διάστατο και τελικά σε 1-διάστατο πρόβλημα. Από τα μονοδιάστατα προφίλ προκύπτουν, με χρήση της θεωρίας μεταβατικών καταστάσεων, οι σταθερές ρυθμού για τις μεταβασεις των ροφούμενων μορίων μεταξύ των θέσεων ρόφησης. Ένα πλεονέκτημα της τροποποιημένης αυτής μεθοδολογίας, που ονομάστηκε Καμπυλόγραμμη Δειγματοληψία Ομπρέλας ή CUS (Curvilinear Umbrella Sampling) σε σχέση με την απλή δειγματοληψία ομπρέλας, είναι ότι μπορεί να εφαρμοστεί και σε μή ευθύγραμμα μονοπάτια μετάβασης. Από τις σταθερές ρυθμού που υπολογιζονται με χρήση της CUS, υπολογίζεται ο συντελεστής αυτoδιάχυσης με χρήση μίας επίσης νέας μεθοδολογίας που εισάγεται σε αυτήν τη διατριβή. Πιο συγκεκριμένα, εφευρέθηκε μια νέα μεθοδολογία για τον υπολογισμό του συντελεστή αυτοδιάχυσης από τις σταθερές ρυθμού των στοιχειωδών αλμάτων και γενικά την ακριβή περιγραφή της μεσοσκοπικής συμπεριφοράς των ροφούμενων μορίων μέσα σε νανοπορώδη υλικά τόσο για μικρούς όσο και για μεγάλους χρόνους. Η μεθοδολογία αυτή είναι εναλλακτική προς την κινητική προσομοίωση Monte Carlo (kinetic Monte Carlo, KMC) και εμφανίζει σαφή πλεονεκτήματα έναντι της τελευταίας. Η νέα μέθοδος επιλύει την εξίσωση εξέλιξης (εξίσωση Master) για ένα σύστημα που προσλαμβάνει διακριτές καταστάσεις που εμφανίζουν περιοδικότητα στο χώρο. Το θεωρούμενο δίκτυο καταστάσεων περιέχει 2ν εικόνες μιας μοναδιαίας κυψελίδας n καταστάσεων, διατεταγμένες κατά μήκος μιας διεύθυνσης του χώρου. Στα άκρα του δικτύου επικρατούν περιοδικές οριακές συνθήκες. Αναλύθηκε η δομή του (2^νn)x(2^νn) συμμετρικού πίνακα των σταθερών ρυθμού και βρέθηκε ένα αναδρομικό σχήμα για τον υπολογισμό των ιδιοτιμών και των ιδιοδιανυσμάτων. Ως αποτέλεσμα εκφράστηκαν αναλυτικά οι πιθανότητες κατάληψης, συναρτήσει του χρόνου, όλων των καταστάσεων που μπορεί να προσλάβει το σύστημα ξεκινώντας από συγκεκριμένες αρχικές συνθήκες. Οι αναλυτικές αυτές εκφρασεις βασίζονται στη διαγωνοποίηση nxn πινάκων οι οποίοι σχηματίζονται από τις καταστάσεις μίας και μόνο μίας κυψελίδας. Αυτή η νέα μέθοδος, εν ονόματι MESoRReD (Master Equation Solution by Rercursive Reduction of Dimensionality) εφαρμόστηκε στο πρόβλημα της διάχυσης το Ξένου στο ζεόλιθο σιλικαλίτη-1 όταν αυτό βρίσκεται σε χαμηλές συγκεντρώσεις και σε χαμηλές θερμοκρασίες, χρησιμοποιώντας τις σταθερες ρυθμού από την δουλειά των June et al., J. Phys. Chem., 95, 8866 (1991). Η μέθοδος MESoRReD μας δίνει ένα τανυστή διαχυτότητας για το συγκεκριμένο σύστημα ο οποίος διαφέρει κατά λιγότερο από 3% από τις τιμές που είχαν παραχθεί στο παρελθόν από προσομοιώσεις Κινητικής Monte Carlo (ΚΜC) αλλά και από προσομοιώσεις ΚΜC που διεξάχθηκαν στα πλαίσια της παρούσας διδακτορική διατριβή. Οι υπολογιστικές απαιτήσεις της συγκεκριμενης μεθόδου συγκρίνονται με εκείνες της ΚΜC, και της αριθμητικής επίλυσης της εξίσωσης Master με την μεθοδο Euler, καθώς και με άμεσες προσομοιώσεις Μοριακής Δυναμικής. Στο πρόβλημα της διάχυσης του Ξένου στο σιλικαλίτη-1, η νέα μέθοδος είναι γρηγορότερη από αυτές τις εναλλακτικές μεθόδους κατά παράγοντες 3.5x10^5, 1.3x10^6, and 1.7x10^7 αντίστοιχα. Ο υπολογιστικός χρόνος που εξοικονομείται από την νέα μέθοδο MESoRReD στη διαγωνοποίηση του πίνακα σταθερών ρυθμού την καθιστά ιδιαίτερα ελκυστική ως μέσο πρόβλεψης δυναμικών φαινομένων μεσαίας αλλά και μεγάλης διάρκειας χρησιμοποιώντας πληροφορίες βασισμένες στο ατομικό επίπεδο. Στην συνέχεια, οι υπολογιστικές μεθοδολογίες που περιγράψαμε εφαρμόστηκαν για την πρόβλεψη του συντελεστή αυτοδιάχυσης του παρα-ξυλολίου και του βενζολίου στο ζεόλιθο σιλικαλίτη-1 όταν τα μόρια αυτά βρίσκονται σε άπειρη αραίωση. Επιπλέον μελετήθηκαν οι κατανομές προσανατολισμού των φαινυλικών δακτυλίων και των μεθυλίων του παραξυλολίου και του βενζολίου ροφημένων στο ζεόλιθο και ελέγχθηκε η ύπαρξη εντροπικών φραγμάτων στη μεταφορική κίνηση. Αναπτύχθηκε μια νέα μέθοδος αναγωγής των καταστάσεων που εμφανίζονται στα προφίλ της τοπικής ελεύθερης ενέργειας και χρησιμοποιήθηκε για τον υπολογισμό των σταθερών ρυθμού στοιχειωδών αλμάτων μέσα στο ζεόλιθο. Τα αποτελέσματα που προέκυψαν για τη διάχυση συμφωνούν με μετρήσεις οιονεί ελαστικής σκέδασης νετρονίων. Το κύριο συμπέρασμα από τις προσομοιώσεις είναι ότι το παραξυλόλιο διαχεέται περίπου 100 φορές γρηγορότερα από το βενζόλιο όταν και τα δύο είναι ροφημένα σε άπειρη αραίωση. Το βενζόλιο υπόκειται σε ισχυρά εντροπικά φράγματα όταν βρίσκεται στην περιοχή διασταύρωσης των καναλιών του σιλικαλίτη, όπου μπορεί να προσλαμβάνει πολλούς προσανατολισμούς. Τα αντίστοιχα φράγματα για το παρα-ξυλόλιο είναι πολύ χαμηλότερα, αντικατοπτρίζοντας τη δυσκολία του κύριου αξονά του να αναπροσανατολιστεί μέσα στις διασταυρώσεις των καναλιών. Εν τέλει εξηγείται για πρώτη φορά το παράδοξο ότι το βαρύτερο και μεγαλύτερο σε μέγεθος μόριο (παρα-ξυλόλιο) διαχέεται γρηγορότερα από ένα ελαφρύτερο και μικρότερο (βενζόλιο) μέσα στο σιλικαλίτη.