Mechanism of ligand exchange studied using transition path sampling

The mechanism of intermolecular ligand exchange has been studied using transition path sampling (TPS) based molecular dynamics (MD) simulations. Specifically, the exchange of solvent molecules bound to unsaturated Cr(CO)5 in methanol solution has been investigated. The results of the TPS simulations have shown that there are multiple steps in the reaction mechanism. The first involves partial dissociation of the coordinated solvent from the Cr metal center followed by association with a new methanol molecule between the normally void first and second solvent layers. After diffusive motion of the exchanging ligands, the last step involves the originally bound methanol molecule moving into the bath continuum followed by solvation of the Cr metal fragment by the exchanging ligand. It has been found that the reaction center (defined as the organometallic fragment and two exchanging ligands only) and the solvent bath have favorable interactions. This is likely due to the adiabatic nature of the ligand exchange transition. The ability to understand the microscopic molecular dynamics of a chemical process based on a free energy analysis is also discussed.     

Efficient transition path sampling for systems with multiple reaction pathways

A new algorithm is developed for sampling transition paths and computing reaction rates. To illustrate the use of this method, we study a two-dimensional system that has two reaction pathways: one pathway is straight with a relatively high barrier and the other is roundabout with a lower barrier. The transition rate and the ratio between the numbers of the straight and roundabout transition paths are computed for a wide range of temperatures. Our study shows that the harmonic approximation for fluctuations about the steepest-descent paths is not valid even at relatively low temperatures and, furthermore, that factors related to entropy have to be determined by the global geometry of the potential-energy surface (rather than just the local curvatures alone) for complex reaction systems. It is reasonable to expect that this algorithm is also applicable to higher dimensional systems.     

Improved transition path sampling methods for simulation of rare events

The free energy surfaces of a wide variety of systems encountered in physics, chemistry, and biology are characterized by the existence of deep minima separated by numerous barriers. One of the central aims of recent research in computational chemistry and physics has been to determine how transitions occur between deep local minima on rugged free energy landscapes, and transition path sampling (TPS) Monte-Carlo methods have emerged as an effective means for numerical investigation of such transitions. Many of the shortcomings of TPS-like approaches generally stem from their high computational demands. Two new algorithms are presented in this work that improve the efficiency of TPS simulations. The first algorithm uses biased shooting moves to render the sampling of reactive trajectories more efficient. The second algorithm is shown to substantially improve the accuracy of the transition state ensemble by introducing a subset of local transition path simulations in the transition state. The system considered in this work consists of a two-dimensional rough energy surface that is representative of numerous systems encountered in applications. When taken together, these algorithms provide gains in efficiency of over two orders of magnitude when compared to traditional TPS simulations.     

Unprejudiced identification of reaction mechanisms from biased transition path sampling

A method for unprejudiced investigation of reaction mechanisms from molecular-dynamics simulations is presented. It combines the transition path sampling approach with a biasing strategy, which (a) allows optimization of transition paths crossing an energy minimum of the transition state surface. The bias is then used to (b) find reaction pathways, which follow different mechanistic routes. In the first step the manifold of similar trajectories that correspond to the same mechanism is reduced to a single characteristic dynamical path. Our method then allows a systematic search for further reaction mechanisms and the related energy barriers. It is illustrated at the example of a single particle in a two-dimensional potential and of the rather complex process of the pressure-induced insertion of a helium atom into a C60 buckyball molecule.     

An efficient transition path sampling algorithm for nanoparticles under pressure

We apply transition path sampling to the simulation of nanoparticles under pressure. As a barostat we use a bath of ideal gas particles that form a stochastically updated atmosphere around the nanoparticle. We justify this algorithm by showing that it preserves the distribution of an ideal gas at constant temperature and pressure by satisfying detailed balance. Based on this result, we present a simple and efficient transition path sampling scheme for the study of activated processes in nanoparticles under pressure. As a first application, we investigate the h-MgO to rocksalt transformation in faceted CdSe nanocrystals. Starting from an artificial mechanism involving a uniform motion of all atoms, trajectories quickly converge towards the dominant mechanism of nucleation and growth along parallel (100) planes.     

Activation energy for a model ferrous-ferric half reaction from transition path sampling

Activation parameters for the model oxidation half reaction of the classical aqueous ferrous ion are compared for different molecular simulation techniques. In particular, activation free energies are obtained from umbrella integration and Marcus theory based thermodynamic integration, which rely on the diabatic gap as the reaction coordinate. The latter method also assumes linear response, and both methods obtain the activation entropy and the activation energy from the temperature dependence of the activation free energy. In contrast, transition path sampling does not require knowledge of the reaction coordinate and directly yields the activation energy [C. Dellago and P. G. Bolhuis, Mol. Simul. 30, 795 (2004)]. Benchmark activation energies from transition path sampling agree within statistical uncertainty with activation energies obtained from standard techniques requiring knowledge of the reaction coordinate. In addition, it is found that the activation energy for this model system is significantly smaller than the activation free energy for the Marcus model, approximately half the value, implying an equally large entropy contribution.     

Exploring reaction pathways with transition path and umbrella sampling: application to methyl maltoside

The transition path sampling (TPS) method is a powerful approach to study chemical reactions or transitional properties on complex potential energy landscapes. One of the main advantages of the method over potential of mean force methods is that reaction rates can be directly accessed without knowledge of the exact reaction coordinate. We have investigated the complementary nature of these two differing approaches, comparing transition path sampling with the weighted histogram analysis method to study a conformational change in a small model system. In this case study, the transition paths for a transition between two rotational conformers of a model disaccharide molecule, methyl beta-D-maltoside, were compared with a free energy surface constrained by the two commonly used glycosidic (phi,psi) torsional angles. The TPS method revealed a reaction channel that was not apparent from the potential of mean force method, and the suitability of phi and psi as reaction coordinates to describe the isomerization in vacuo was confirmed by examination of the transition path ensemble. Using both transition state theory and transition path sampling methods, the transition rate was estimated. We have estimated a characteristic time between transitions of approximately 160 ns for this rare isomerization event between the two conformations of the carbohydrate. We conclude that transition path sampling can extract subtle information about the dynamics not apparent from the potential of mean force method. However, in calculating the reaction rate, the transition path sampling method required 27.5 times the computational effort than was needed by the potential of mean force method.     

A transition path sampling study of the reaction catalyzed by the enzyme chorismate mutase

The study of the chemical steps in enzyme-catalyzed reactions represents a challenge for molecular simulation techniques. One concern is how to calculate paths for the reaction. Common techniques include the definition of a reaction coordinate in terms of a small set of (normally) geometrical variables or the determination of minimum energy paths on the potential energy surface of the reacting system. Both have disadvantages, the former because it presupposes knowledge of which variables are likely to be important for reaction and the latter because it provides a static picture and dynamical effects are ignored. In this paper, we employ the transition path sampling method developed by Chandler and co-workers, which overcomes some of these limitations. The reaction that we have chosen is the chorismate-mutase-catalyzed conversion of chorismate into prephenate, which has become something of a test case for simulation studies of enzyme mechanisms. We generated an ensemble of approximately 1000 independent transition paths for the reaction in the enzyme and another approximately 500 for the corresponding reaction in solution. A large variety of analyses of these paths was performed, but we have concentrated on characterizing the transition state ensemble, particularly the flexibility of its structures with respect to other ligands of the enzyme and the time evolution of various geometrical and energetic properties as the reaction proceeds. We have also devised an approximate technique for locating transition state structures along the paths.     

How enzyme dynamics helps catalyze a reaction in atomic detail: a transition path sampling study

We have applied the Transition Path Sampling algorithm to the reaction catalyzed by the enzyme Lactate Dehydrogenase. This study demonstrates the ease of scaling Transition Path Sampling for applications on many degree of freedom systems, whose energy surface is a complex terrain of valleys and saddle points. As a Monte Carlo importance sampling method, transition path sampling is capable of surmounting barriers in path phase space and focuses simulation on the rare event of enzyme catalyzed atom transfers. Generation of the transition path ensemble, for this reaction, resolves a paradox in the literature in which some studies exposed the catalytic mechanism of hydride and proton transfer by lactate dehydrogenase to be concerted and others stepwise. Transition path sampling has confirmed both mechanisms as possible paths from reactants to products. With the objective to identify a generalized, reduced reaction coordinate, time series of both donor-acceptor distances and residue distances from the active site have been examined. During the transition from pyruvate to lactate, residues located behind the transferring hydride collectively compress toward the active site causing residues located behind the hydride acceptor to relax away. It is demonstrated that an incomplete compression/relaxation transition across the donor-acceptor axis compromises the reaction.     

Generation of initial trajectories for transition path sampling of chemical reactions with ab initio molecular dynamics

Transition path sampling is an innovative method for focusing a molecular dynamics simulation on a reactive event. Although transition path sampling methods can generate an ensemble of reactive trajectories, an initial reactive trajectory must be generated by some other means. In this paper, the authors have evaluated three methods for generating initial reactive trajectories for transition path sampling with ab initio molecular dynamics. The authors have tested each of these methods on a set of chemical reactions involving the breaking and making of covalent bonds: the 1,2-hydrogen elimination in the borane-ammonia adduct, a tautomerization, and the Claisen rearrangement. The first method is to initiate trajectories from the potential energy transition state, which was effective for all reactions in the test set. Assigning atomic velocities found using normal mode analysis greatly improved the success of this method. The second method uses a high temperature molecular dynamics simulation and then iteratively reduces the total energy of the simulation until a low temperature reactive trajectory is found. This was effective in generating a low temperature trajectory from an initial trajectory run at 3000 K of the tautomerization reaction, although it failed for the other two. The third uses an orbital based bias potential to find a reactive trajectory and uses this trajectory to initiate an unbiased trajectory. The authors found that a highest occupied molecular orbital-lowest unoccupied molecular orbital bias could be used to find a reactive trajectory for the Claisen rearrangement, although it failed for the other two reactions. These techniques will help make it practical to use transition path sampling to study chemical reaction mechanisms that involve bond breaking and forming.     

Optimizing transition interface sampling simulations

We demonstrate that a recently proposed adaptive optimization algorithm for forward flux sampling simulations [E. E. Borrero and F. A. Escobedo, J. Chem. Phys. 129, 024115 (2008)] can be easily applied within the framework of transition interface sampling. This optimization algorithm systematically identifies the kinetic bottlenecks along the order parameter used to partition phase space via interfaces and improves the statistical accuracy of the reaction rate constant estimate. In different versions of the algorithm, the number or the placement of the interfaces (or both) are varied in order to allocate the numerical effort in a balanced way. The algorithm is demonstrated for a simple two-dimensional model and for the dipole flip transition of icelike structures inside carbon nanotubes. For these test systems, the optimization yielded an efficiency increase by a factor of 2-15.     

Selective sampling of transition paths

In this short paper, we introduce an approximate method for the quick estimate of rate constants based on a simple sampling method of reactive transition paths over high energy barriers. It makes use of the previously introduced accelerated molecular dynamics (MD) simulation method to generate initial points for trajectory shooting. The accelerated MD simulations, although with the loss of real dynamics, lead to a quick calculation of thermodynamic properties and at the same time produce an ensemble of configurations with an enhanced sampling over the phase space that is more "reactive." The forward/backward trajectory shooting as that used in the transition path sampling method is then initiated from the configurations obtained from accelerated MD simulations to generate transition paths on the original unbiased potential. This method selectively enhances sampling of successful trajectories and at the same time accelerates significantly the calculation of rate constants.     

Unbiased sampling of lattice Hamilton path ensembles

Hamilton paths, or Hamiltonian paths, are walks on a lattice which visit each site exactly once. They have been proposed as models of globular proteins and of compact polymers. A previously published algorithm [Mansfield, Macromolecules 27, 5924 (1994)] for sampling Hamilton paths on simple square and simple cubic lattices is tested for bias and for efficiency. Because the algorithm is a Metropolis Monte Carlo technique obviously satisfying detailed balance, we need only demonstrate ergodicity to ensure unbiased sampling. Two different tests for ergodicity (exact enumeration on small lattices, nonexhaustive enumeration on larger lattices) demonstrate ergodicity unequivocally for small lattices and provide strong support for ergodicity on larger lattices. Two other sampling algorithms [Ramakrishnan et al., J. Chem. Phys. 103, 7592 (1995); Lua et al., Polymer 45, 717 (2004)] are both known to produce biases on both 2x2x2 and 3x3x3 lattices, but it is shown here that the current algorithm gives unbiased sampling on these same lattices. Successive Hamilton paths are strongly correlated, so that many iterations are required between statistically independent samples. Rules for estimating the number of iterations needed to dissipate these correlations are given. However, the iteration time is so fast that the efficiency is still very good except on extremely large lattices. For example, even on lattices of total size 10x10x10 we are able to generate tens of thousands of uncorrelated Hamilton paths per hour of CPU time.     

Overcoming barriers in trajectory space: mechanism and kinetics of rare events via Wang-Landau enhanced transition path sampling

Within the framework of transition path sampling (TPS), activation energies can be computed as path ensemble averages without a priori information about the reaction mechanism [C. Dellago and P. G. Bolhuis, Mol. Simul. 30, 795 (2004)]. Activation energies computed for different conditions can then be used to determine by numerical integration the rate constant for a system of interest from the rate constant known for a reference system. However, in systems with complex potential energy surfaces, multiple reaction pathways may exist making ergodic sampling of trajectory space difficult. Here, we present a combination of TPS with the Wang-Landau (WL) flat-histogram algorithm for an efficient sampling of the transition path ensemble. This method, denoted by WL-TPS, has the advantage that from one single simulation, activation energies at different temperatures can be determined even for systems with multiple reaction mechanisms. The proposed methodology for rate constant calculations does not require the knowledge of the reaction coordinate and is generally applicable to Arrhenius and non-Arrhenius processes. We illustrate the applicability of this technique by studying a two-dimensional toy system consisting of a triatomic molecule immersed in a fluid of repulsive soft disks. We also provide an expression for the calculation of activation volumes from path averages such that the pressure dependence of the rate constant can be obtained by numerical integration.     

Transition path sampling study of classical rate-promoting vibrations

It is now widely accepted that there is a class of enzymatic proton transfer reactions, which proceed through quantum tunneling. In a series of papers we have argued that some experimental features of these reactions can be explained by assuming the presence of a "rate-promoting" vibration which brings donor and acceptor closer together, thus leading to rate enhancement. There has never been a study of this effect for classical systems. We used transition path sampling to study the equivalent classical problem and found a complicated dynamical behavior that cannot be captured by transition state theory. Slow promoting vibrations lead to reactive trajectories that overshoot the saddle point, but on the other hand the short period of fast oscillations allows the reactants to stay only briefly in a low-barrier regime. There is a competition between these effects, which results to an intermediate value for the frequency of the rate-promoting vibration that is optimal for enhancing the rate.     

The dynamics of water exchange in gadolinium DOTA complexes studied by transition path sampling and potential of mean force methods

The mechanism of water exchange at the Gd centre of the two isomers of [Gd(iii)DOTA](-) (gadolinate(1-), [1,4,7,10-tetraazacyclododecane-1,4,7,10-tetraacetato(4-)-N1,N4,N7,N10,O1,O4,O7,O10]) has been explored using transition path sampling and potential of mean force methods to sample those regions of phase space inaccessible to standard molecular dynamics simulation. We find that there are definite differences in the details of the solvent rearrangement accompanying the exchange of the capping water molecule for the two isomers. We conclude that these solvent effects, rather than any differences in the binding energy of the capping water, are central in determining the exchange rate. We find that the potential of mean force studies yield absolute and relative rates of water exchange for the two isomers that are in good agreement with experiment.     

Path ensembles and path sampling in nonequilibrium stochastic systems

Markovian models based on the stochastic master equation are often encountered in single molecule dynamics, reaction networks, and nonequilibrium problems in chemistry, physics, and biology. An efficient and convenient method to simulate these systems is the kinetic Monte Carlo algorithm which generates continuous-time stochastic trajectories. We discuss an alternative simulation method based on sampling of stochastic paths. Utilizing known probabilities of stochastic paths, it is possible to apply Metropolis Monte Carlo in path space to generate a desired ensemble of stochastic paths. The method is a generalization of the path sampling idea to stochastic dynamics, and is especially suited for the analysis of rare paths which are not often produced in the standard kinetic Monte Carlo procedure. Two generic examples are presented to illustrate the methodology.     

Rate constants for diffusive processes by partial path sampling

We introduce a path sampling method for the computation of rate constants for complex systems with a highly diffusive character. Based on the recently developed transition interface sampling (TIS) algorithm this procedure increases the efficiency by sampling only parts of complete transition trajectories. The algorithm assumes the loss of memory for diffusive progression along the reaction coordinate. We compare the new partial path technique to the TIS method for a simple diatomic system and show that the computational effort of the new method scales linearly, instead of quadratically, with the width of the diffusive barrier. The validity of the memory loss assumption is also discussed.     

Variational transition state theory without the minimum-energy path

In this paper we propose a method for carrying out variational transition state theory calculations without first obtaining a converged minimum-energy path (MEP). We illustrate the method in two ways, first of all by employing an unconverged MEP and secondly by using a dynamically optimized distinguished reaction path. Preliminary tests of the algorithm for the reactions OH+H₂→H₂O+H and C₂H₅→C₂H₄+H are very encouraging. Received: 22 January 1997 / Accepted: 11 March 1997     

Single path transition probabilities in potential-curve crossing phenomena

For a potential curve crossing with radial coupling (Landau-Zener crossing) we have investigated the transition probabilities corresponding to a single passage through the crossing point. Special emphasis has been given to the convergence behaviour as a function of the width of the coupling region. It has been found that the convergence is much faster, i.e., the effective width is smaller, for the adiabatic representation of the system than for the diabatic one. This is true even for very high energies when the system has a large probability to remain in the same diabatic state during its passage through the crossing point. For the example of the curve crossing at 7 Å for Na I surprisingly good agreement is found between the Landau-Zener formula and our exact numerical result.