Developing the molecular modelling of diffusion in zeolites as a high throughput catalyst screening technique

Molecular modelling techniques have been used to screen zeolite catalysts for their suitability for organic synthesis. For example, we have used these techniques to study the alkylation of aromatic molecules which are important in the fine-chemical and drug industries. A survey of all such efforts is reviewed in this article. The application of molecular modelling techniques in a systematic manner is an efficient first step in the design of zeolite catalysts. As a qualitative screening tool, molecular graphics is used to visualize how well the reactant and product molecules fit inside the pores of the zeolites. Using a hybrid of several molecular modelling methods, which combines molecular dynamics (MD) and Monte Carlo methods with energy minimization, it is possible to determine the minimum energy locations of the molecules inside the zeolites cages. The minimum energy configurations determined by this hybrid method are taken as a starting point for diffusion of the molecules through the zeolite channels. When a molecule is allowed to diffuse through zeolite channel, the molecule attains some maxima and minima in its diffusion energy profile. From the differences between a maximum and a minimum energy configuration, the diffusion energy barrier for the molecule can be calculated in the zeolites. By comparing the diffusion energy barriers for various isomers of a molecule in different zeolites, it is possible to find out the most suitable zeolite for achieving the required shape-selectivity. In addition, factors influencing the diffusivity of the molecules and consequently the shape selectivity are derived. The list of factors and their relative importance are analysed to derive valuable guidelines to design shape-selective zeolite catalysts for a given reaction. Thus, the ultimate aim of these studies is to develop a high throughput computational screening process for the selection of shape-selective zeolite catalysts for various reactions. The dynamic behaviour of molecules inside the pores of zeolites can be studied using MD methods. Since MD is computationally time consuming, it is more efficient to screen the possible zeolite catalysts by energy minimization methods and then perform MD in specific zeolites. More accurate values of diffusivity of the molecules can be calculated using MD methods, and these values can be correlated with the shape-selectivity observed experimentally and /or derived from diffusion energy barrier calculations.     

Computation of binding free energy with molecular dynamics and grand canonical Monte Carlo simulations

The binding of a ligand to a receptor is often associated with the displacement of a number of bound water molecules. When the binding site is exposed to the bulk region, this process may be sampled adequately by standard unbiased molecular dynamics trajectories. However, when the binding site is deeply buried and the exchange of water molecules with the bulk region may be difficult to sample, the convergence and accuracy in free energy perturbation (FEP) calculations can be severely compromised. These problems are further compounded when a reduced system including only the region surrounding the binding site is simulated. To address these issues, we couple molecular dynamics (MD) with grand canonical Monte Carlo (GCMC) simulations to allow the number of water to fluctuate during an alchemical FEP calculation. The atoms in a spherical inner region around the binding pocket are treated explicitly while the influence of the outer region is approximated using the generalized solvent boundary potential (GSBP). At each step during thermodynamic integration, the number of water in the inner region is equilibrated with GCMC and energy data generated with MD is collected. Free energy calculations on camphor binding to a deeply buried pocket in cytochrome P450cam, which causes about seven water molecules to be expelled, are used to test the method. It concluded that solvation free energy calculations with the GCMC/MD method can greatly improve the accuracy of the computed binding free energy compared to simulations with fixed number of water.     

Theoretical modeling of enzyme reactions: the thermodynamics of formation of compound 0 in horseradish peroxidase

In this paper, by using the perturbed matrix method (PMM) in combination with basic statistical mechanical relations both based on nanosecond time-scale molecular dynamics (MD) simulations, we quantitatively address the thermodynamics of compound 0 (Cpd 0) formation in horseradish peroxidase (HRP) enzyme. Our results, in the same trend of low-temperature experimental data, obtained in cryoenzymology studies indicate that such a reaction can be described essentially as a stepwise spontaneous process: a first step mechanically constrained, strongly exothermic proton transfer from the heme-H2O2 complex to the conserved His42, followed by a solvent-protein relaxation involving a large entropy increase. Critical evaluation of PMM/MD data also reveals the crucial role played by specific residues in the reaction pocket and, more in general, by the conformational fluctuations of the overall environment in physiological conditions.     

A theoretical analysis of rate constants and kinetic isotope effects corresponding to different reactant valleys in lactate dehydrogenase

In some enzymatic systems large conformational changes are coupled to the chemical step, in such a way that dispersion of rate constants can be observed in single-molecule experiments, each corresponding to the reaction from a different reactant valley. Under this perspective here we present a computational study of pyruvate to lactate transformation catalyzed by lactate dehydrogenase. The reaction consists of a hydride transfer and a proton transfer that seem to take place concertedly. The degree of asynchronicity and the energy barrier depend on the particular starting reactant valley. In order to estimate rate constants we used a free energy perturbation technique adapted to follow the intrinsic reaction coordinate for several possible reaction paths. Tunneling effects are also obtained with a slightly modified version of the ensemble-averaged variational transition state theory with multidimensional tunneling contributions. According to our results the closure of the active site by means of a flexible loop can lead to the formation of different reactant complexes displaying different features in the disposition of some key residues (such as Arg109), interactions with the substrate and number of water molecules in the active site. The chemical step of the reaction takes place with a different reaction rate from each of these complexes. Finally, primary kinetic isotope effects for replacement of the transferring hydrogen of the cofactor with a deuteride are in good agreement with experimental observations, thus validating our methodology.     

Quantum chemical prediction of reaction pathways and rate constants for dissociative adsorption of CO(x) and NO(x) on the graphite (0001) surface

We present predictions of reaction rate constants for dissociative adsorption reactions of CO(x) (x = 1, 2) and NO(x) (x = 1, 2) molecules on the basal graphite (0001) surface based on potential energy surfaces (PES) obtained by the integrated ONIOM(B3LYP:DFTB-D) quantum chemical hybrid approach with dispersion-augmented density functional tight binding (DFTB-D) as low level method. Following an a priori methodology developed in a previous investigation of water dissociative adsorption reactions on graphite, we used a C(94)H(24) dicircumcoronene graphene slab as model system for the graphite surface in finite-size molecular structure investigations, and single adsorbate molecules reacting with the pristine graphene sheet. By employing the ONIOM PES information in RRKM theory we predict reaction rate constants in the temperature range between 1,000 and 5,000 K. We find that among CO(x) and NO(x) adsorbate species, the dissociative adsorption reactions of CO(2) and both radical species NO and NO(2) are likely candidates as a cause for high temperature oxidation and erosion of graphite (0001) surfaces, whereas reaction with CO is not likely to lead to long-lived surface defects. High temperature quantum chemical molecular dynamics simulations (QM/MD) at T = 5,000 K using on-the-fly DFTB-D energies and gradients confirm the results of our PES study.     

The generalized transition state method

Current theories of unimolecular reaction rates are based on the transition state method which replaces internal reactant dynamics by an assumption of internal equilibrium. The present work is devoted to the development of generalized transition state method which allows effects such as nonergodicity and non-exponential decay to be accounted for within a simple theoretical framework. The derivation is quantum mechanical and not limited by any weak perturbation assumption. An effective hamiltonian is constructed for the reactant dynamics. The loss of amplitude due to reaction is accounted for by a dissipative term in the hamiltonian which is obtained on a phenomenological basis. The diagonalization of the hamiltonian allows the decay of reactant state to be predicted. The decay information is then used to set up a non-markovian master equation which in turn yields the rate coefficient for the reaction. The accuracy of the method is tested in one-dimensional model calculations in which particular attention is paid to decay by quantum mechanical tunneling through a potential barrier.     

Cavity Catalysis by Cooperative Vibrational Strong Coupling of Reactant and Solvent Molecules

Here, we report the catalytic effect of vibrational strong coupling (VSC) on the solvolysis of para‐nitrophenyl acetate (PNPA), which increases the reaction rate by an order of magnitude. This is observed when the microfluidic Fabry–Perot cavity in which the VSC is generated is tuned to the C=O vibrational stretching mode of both the reactant and solvent molecules. Thermodynamic experiments confirm the catalytic nature of VSC in the system. The change in the reaction rate follows an exponential relation with respect to the coupling strength of the solvent, indicating a cooperative effect between the solvent molecules and the reactant. Furthermore, the study of the solvent kinetic isotope effect clearly shows that the vibrational overlap of the C=O vibrational bands of the reactant and the strongly coupled solvent molecules is critical for the catalysis in this reaction. The combination of cooperative effects and cavity catalysis confirms the potential of VSC as a new frontier in chemistry.     

Molecular beam measurements of small specific rate constants for unimolecular reactions

The investigation of unimolecular reactions with small rate constants is difficult owing to competing processes (inelastic collisions and bimolecular reactions) and the diffusion of reactant and product molecules out of the detection volume. For this reason, a new experimental approach for the measurement of specific rate constants in a molecular beam experiment has been exploited; instead of monitoring the temporal change of intensity as in a cell experiment, we monitor the spatial change along the molecular beam axis after laser excitation. For a given particle velocity the flight path between excitation and detection region defines the reaction time. By varying the distance the specific rate constant can be determined directly both from the decrease in the number density of reactant molecules as well as from the increase in product molecules. As a model system, the laser-induced (λ = 193 nm) photodissociation of mesitylene (trimethylbenzene) is studied. Previous experiments on the specific rate constant of mesitylene at this excitation energy differ between each other by about a factor of ten. By combining the new results with measurements at higher excitation energies, rate constants over a range of two orders of magnitude are now available for this reaction. The differences between the various experimental results are discussed within the framework of a statistical theory.     

ReaxFF MD局部区域反应追踪与物理性质可视化分析

Recently, the application of ReaxFF based reactive molecular dynamics simulation (ReaxFF MD) in complex processes of pyrolysis, oxidation and catalysis has attracted considerable attention. The analysis of the simulation results of these processes is challenging owing to the complex chemical reactions involved, coupled with their dynamic physical properties. VARxMD is a leading tool for the chemical reaction analysis and visualization of ReaxFF MD simulations, which allows the automated analysis of reaction sites to get overall reaction lists, evolution trends of reactants and products, and reaction networks of specified reactants and products. The visualization of the reaction details and dynamic evolution profiles are readily available for each reactant and product. Additionally, the detailed reaction sites of bond breaking and formation are available in 2D chemical structure diagrams and 3D structure views; for specified reactions, they are categorized on the basis of the chemical structures of the bonding sites or function groups in the reacting species. However, the current VARxMD code mainly focuses on global chemical reaction information in the simulation system of the ReaxFF MD, and is incapable of locally tracking the chemical reaction and physical properties in a 3D picked zone. This work extends the VARxMD from global analysis to a focused 3D zone picked interactively from the 3D visualization modules of VARxMD, as well as physical property analysis to complement reaction analysis. The analysis of reactions and physical properties can be implemented in three steps: picking and drawing a 3D zone, identifying molecules in the picked zone, and analyzing the reactions and physical properties of the picked molecules. A 3D zone can be picked by specifying the geometric parameters or drawing on a screen using a mouse. The picking of a cuboid or sphere was implemented using the VTK 3D view libraries by specifying geometric parameters. The interactive 3D zone picking was implemented using a combination of observer and command patterns in the VTK visualization paradigm. The chemical reaction tracking and dynamic radial distribution function (RDF) of the 3D picked zone was efficiently implemented by inheriting data obtained from the global analysis of VARxMD. The reaction tracking between coal particles in coal pyrolysis simulation and dynamic structure characterization of carbon rich cluster formation in the thermal decomposition of an energetic material are presented as application examples. The obtained detailed reactions between the coal particles and comparison of the reaction between the locally and globally picked areas in the cuboid are helpful in understanding the role of micro pores in coal particles. The carbon to carbon RDF analysis and comparison of the spherical region picked for the layered molecular clusters in the pyrolysis system of the TNT crystal model with the standard RDF of the 5-layer graphene demonstrate the extended VARxMD as a chemical structure characteristic tool for detecting the dynamic formation profile of carbon rich clusters in the pyrolysis of TNT. The extended capability of VARxMD for a 3D picked zone of a ReaxFF MD simulation system can be useful for interfacial reaction analysis in a catalysis system, hot spot formation analysis in the detonation of energetic material systems, and particularly the pyrolysis or oxidation processes of coal, biomass, polymers, hydrocarbon fuels, and energetic materials.     

The midpoint method for parallelization of particle simulations

The evaluation of interactions between nearby particles constitutes the majority of the computational workload involved in classical molecular dynamics (MD) simulations. In this paper, we introduce a new method for the parallelization of range-limited particle interactions that proves particularly suitable to MD applications. Because it applies not only to pairwise interactions but also to interactions involving three or more particles, the method can be used for evaluation of both nonbonded and bonded forces in a MD simulation. It requires less interprocessor data transfer than traditional spatial decomposition methods at all but the lowest levels of parallelism. It gains an additional practical advantage in certain commonly used interprocessor communication networks by distributing the communication burden more evenly across network links and by decreasing the associated latency. When used to parallelize MD, it further reduces communication requirements by allowing the computations associated with short-range nonbonded interactions, long-range electrostatics, bonded interactions, and particle migration to use much of the same communicated data. We also introduce certain variants of this method that can significantly improve the balance of computational load across processors.     

Ab initio MD simulations of a prototype of methyl chloride hydrolysis with explicit consideration of three water molecules: a comparison of MD trajectories with the IRC path

Ab initio molecular dynamics simulations at the Hartree-Fock/6-31G level of theory are performed on methyl chloride hydrolysis with explicit consideration of one solute and two solvent water molecules at a temperature of 298 K. The reaction involves the formation of a reactant complex and the energy surface to the transition state is found to be simple. Two types of trajectories toward the product are observed. In the first type, the system reaches an intermediate complex (complex-P1) region after two nearly concerted proton transfers involving the attacking water molecule and the solvent water molecules. These trajectories resemble the intrinsic reaction coordinate trajectory. The thermal motion of the atoms leads the system to another intermediate complex (complex-P2) region. A second type of trajectory is found in which the system reaches the complex-P2 region directly after the proton transfers. In both of these forward trajectories, back proton transfers lead the system to a final complex-F region which resembles protonated methanol. Received: 3 July 1998 / Accepted: 2 September 1998 / Published online: 15 February 1999     

Time series analysis of ion dynamics in glassy ionic conductors obtained by a molecular dynamics simulation

We present several characteristics of ionic motion in glassy ionic conductors brought out by time series analysis of molecular dynamics (MD) simulation data. Time series analysis of data obtained by MD simulation can provide crucial information to describe, understand and predict the dynamics in many systems. The data have been treated by the singular spectrum analysis (SSA), which is a method to extract information from noisy short time series and thus provide insight into the unknown or partially unknown dynamics of the underlying system that generated the time series. Phase-space plot reconstructed using the principal components of SSA exhibited complex but clear structures, suggesting the deterministic nature of the dynamics.     

Structural dissimilarity sampling with dynamically self‐guiding selection

Structural dissimilarity sampling (SDS) has been proposed as an enhanced conformational sampling method for reproducing the structural transitions of a given protein. SDS consists of cycles of two steps: (1) Selections of initial structures with structural dissimilarities by referring to a measure. (2) Conformational resampling by restarting short‐time molecular dynamics (MD) simulations from the initial structures. In the present study, an efficient measure is proposed as a dynamically self‐guiding selection to accelerate the structural transitions from a reactant state to a product state as an extension to the original SDS. In the extended SDS, the inner product (IP ) between the reactant and the snapshots generated by short‐time MD simulations are evaluated and ranked according to the IP s at every cycle. Then, the snapshots with low IP s are selected as initial structures for the short‐time MD simulations. This scheme enables one to choose dissimilar and distant initial structures from the reactant, and thus the initial structures dynamically head towards the product, promoting structural transitions from the reactant. To confirm the conformational sampling efficiency, the extended SDS was applied to maltodextrin binding protein (MBP), and we successfully reproduced the structural transition from the open to closed states with submicrosecond‐order simulation times. However, a conventional long‐time MD simulation failed to reproduce the same structural transition. We also compared the performance with that obtained by the ordinary SDS and other sampling techniques that have been developed by us to characterize the possible utility of the extended SDS for actual applications. © 2017 Wiley Periodicals, Inc.     

A quantum theory of chemical processes and reaction rates based on diabatic electronic functions coupled in an external field

The method discussed in this work provides a theoretical framework where simple chemical reactions resemble any other standard quantum process, i.e., a transition in quantum state mediated by the electromagnetic field. In our approach, quantum states are represented as a superposition of electronic diabatic basis functions, whose amplitudes can be modulated by the field and by the external control of nuclear configurations. Using a one-dimensional three-state model system, we show how chemical structure and dynamics can be represented in terms of these control parameters, and propose an algorithm to compute the reaction probabilities. Our analysis of effective energy barriers generalizes previous ideas on structural similarity between reactant, and product, and transition states using the geometry of conventional reaction paths. In the present context, exceptions to empirical rules such as the Hammond postulate appear as effects induced by the environment that supplies the external field acting on the quantum system.     

A surface hopping method for chemical reaction dynamics in solution described by diabatic representation: an analysis of tunneling and thermal activation

We present a surface hopping method for chemical reaction in solution based on diabatic representation, where quantum mechanical time evolution of the vibrational state of the reacting nuclei as well as the reaction-related electronic state of the system are traced simultaneously together with the classical motion of the solvent. The method is effective in describing the system where decoherence between reactant and product states is rapid. The diabatic representation can also give a clear picture for the reaction mechanism, e.g., thermal activation mechanism and a tunneling one. An idea of molecular orbital theory has been applied to evaluate the solvent contribution to the electronic coupling which determines the rate of reactive transition between the reactant and product potential surfaces. We applied the method to a model system which can describe complex chemical reaction of the real system. Two numerical examples are presented in order to demonstrate the applicability of the present method, where the first example traces a chemical reaction proceeded by thermal activation mechanism and the second examines tunneling mechanism mimicking a proton transfer reaction.     

A combined ONIOM quantum chemical-molecular dynamics study of zinc-uracil bond breaking in yeast cytosine deaminase

A QM/MM method that combines ONIOM quantum chemistry and molecular dynamics is developed and applied to a step in the deamination of cytosine to uracil in yeast cytosine deaminase (yCD). A two-layer ONIOM calculation is used for the reaction complex, with an inner part treated at a high level for the chemical reaction (bond breaking) and a middle part treated at a lower level for relevant protein residues that are frozen in the quantum optimization. An outer layer (protein and solvent) is treated using MD. Configurations for the entire system are generated using MD and optimized with ONIOM. The method permits the use of high-level quantum calculations along with sufficient configurational sampling to approximate the potential of mean force for certain bond-breaking reactions. A previously proposed reaction mechanism for deamination (Sklenak, S.; Yao, L. S.; Cukier, R. I.; Yan, H. G. J. Am. Chem. Soc. 2004, 126, 14879) requires breaking the bond between a catalytic zinc and the O4 of uracil in order to permit product release. Using an ONIOM approach, direct bond cleavage was found to be energetically unfavorable. In the work presented here, the combined ONIOM MD method is used to show that the barrier for bond cleavage is small, approximately 3 kcal/mol, and, consequently, should not be the rate-limiting step in the reaction.     

Why and how do systems react in thermally fluctuating environments?

Many chemical reactions, including those of biological importance, take place in thermally fluctuating environments. Compared to isolated systems, there arise markedly different features due to the effects of energy dissipation through friction and stochastic driving by random forces reflecting the fluctuation of the environment. Investigation of how robustly the system reacts under the influence of thermal fluctuation, and elucidating the role of thermal fluctuation in the reaction are significant subjects in the study of chemical reactions. In this article, we start with overviewing the generalized Langevin equation (GLE), which has long been used and continues to be a powerful tool to describe a system surrounded by a thermal environment. It has been also generalized further to treat a nonstationary environment, in which the conventional fluctuation-dissipation theorem no longer holds. Then, within the framework of the Langevin equation we present a method recently developed to extract a new reaction coordinate that is decoupled from all the other coordinates in the region of a rank-one saddle linking the reactant and the product. The reaction coordinate is buried in nonlinear couplings among the original coordinates under the influence of stochastic random force. It was ensured that the sign of this new reaction coordinate (= a nonlinear functional of the original coordinates, velocities, friction, and random force) at any instant is sufficient to determine in which region, the reactant or the product, the system finally arrives. We also discuss how one can extend the method to extract such a coordinate from the GLE framework in stationary and nonstationary environments, where memory effects exist in dynamics of the reaction.     

A versatile AMBER-Gaussian QM/MM interface through PUPIL

The PUPIL package (Program for User Package Interfacing and Linking) originally was developed to interface different programs for multiscale calculations in materials sciences (Torras et al., J Comput Aided Mater Des 2006, 13, 201; Torras et al., Comput Phys Commun 2007, 177, 265). Here we present an extension of PUPIL to computational chemistry by interfacing two widely used computational chemistry programs: AMBER (molecular dynamics) and Gaussian (quantum chemistry). The benefit is to allow the application of the advanced MD techniques available in AMBER to a hybrid QM/MM system in which the forces and energy on the QM part can be computed by any of the methods available in Gaussian. To illustrate, we present two example applications: A MD calculation of alanine dipeptide in explicit water, and a use of the steered MD capabilities in AMBER to calculate the free energy of reaction for the dissociation of Angeli's salt.