Finding Transition Pathways in Complex Systems: Throwing Ropes Over Rough Mountain Passes, In The Dark

During the last few years, our research group has developed a general computational method for finding the transition pathways for infrequent events. The algorithm requires no preconceived notion of mechanism or transition state. Called "transition path sampling," the method is metaphorically akin to throwing ropes over rough mountain passes, in the dark. "Throwing ropes" in the sense that one shoots short trajectories, attempting to reach one stable state from another. "In the dark" because high dimensional systems are so complex, it is generally impossible to literally visualize the topography of relevant energy surfaces. In such cases, it is unlikely that the first throw of the rope will be successful, but one can learn from failures, and there should be an optimum algorithm, i.e., sequence of throws, with which success is obtained efficiently. We have discovered and demonstrated this type of sequence, and we believe it opens the way for many heretofore impossible computational studies of the dynamic pathways of chemical reactions in clusters and in condensed phases.

 

Fig. 1: Schematical depiction of the the potential energy of a complex system. Even though such an energy landscape is dense in saddle points only a few of them are relevant for transitions between different basins of attractions. At finite temperature all details of the surface smaller than kT are of minor importance. Since the transition path sampling method does not rely on identifying saddle points in the potential energy surface, it is the tool of choice to study transitions in complex systems.

Rare but important events in complex systems

Often, dynamical processes of interest occur on time scales that are very long compared to the fundamental unit of time. For example, the dissociation of a weak acid in water might occur with a half life of, say, a milli second (i.e., 10^-3 s), while elementary steps of molecular motions in water occur in femto seconds (i.e., 10^-15 s). Similarly, time scales for folding for the smallest of proteins is in the range of milli to micro seconds, while that for small amplitude motions of amino acid side chains and water solvent is again femto seconds.

This wide disparity of time scales can present serious computational challenges. For instance, consider a computed trajectory for a system containing a weak acid molecule and a bath of a few hundred water molecules. Within one or two orders of magnitude -- depending on computing equipment and algorithm -- 1 s of computation time is required to advance the system for what would correspond to 1 fs of physical time. Therefore, typically 10¹² s of computing time seems to be required to find one example of an event leading to acid dissociation. A representative sampling of pathways to dissociation would therefore seem to be an impossibly impractical computational task.

One way to get around this problem is to focus on the dynamical bottleneck for the rare event -- the transition state. If the location of its threshold is known, one may construct a scheme where the system is first moved reversibly to that state, and then many fleeting trajectories can be initiated from that state. The first step determines the reversible work and thus the probability for being at the transition state, and the subsequent trajectories determine the probability for successfully crossing the threshold. Together, they give the rate for the rare event [1]. While theoretically sound, this two step procedure is highly limited in applicability because it presupposes knowledge of the transition state. In most interesting cases, transition states are not known.

For low dimensional systems involving, say three or fewer atoms, transition states can be located numerically with various algorithms that systematically search for saddle points in the potential energy surface. For higher dimensional systems, however, the potential energy surface will typically contain many saddle points. Explicit enumeration of saddle points is feasible for a cluster of ten or fewer atoms, but provides no means to distinguish saddle points that are dynamically irrelevant from those that are dynamically relevant. For a chaotic system -- large polyatomic molecules, large clusters and condensed phases -- potential energy surfaces are rough on the scale of thermal energies, kT, and dense in saddle points. Effectively, therefore, there is generally an uncountable number of transition states. Searching for a few such states is insufficient. Instead, one wants to locate and sample an ensemble of transition states. We have discovered how this sampling can be accomplished.

Transition path sampling

The basic idea is a generalization of standard Monte Carlo procedures. In its standard form, a Monte Carlo calculation performs a random walk in configuration space. The walk is biased to ensure that the most important regions of configuration space are adequately sampled. Specifically, in a Monte Carlo random walk, configuration x is visited in proportion to its probability p(x). The walk may be initiated far from a typical configuration, but after some equilibration period, the bias drives the system to the important regions of configuration space. This feature, crucial to its success, is called "importance sampling."

 

Fig. 2: In a Mon te-Carlo simulation one generates a random walk in configuration space according to the probablity distribution p(x)=exp(-V(x)/kT). Along this walk a new configuration x' is generated by displacing the old configuration x a little. Then x' is accepted or rejected. If the step goes downhill in energy, i.e. if the new configuration has a higher probability than the old one, x' is always accepted. Uphill moves, on the other hand, are only accepted with a probability p(x')/p(x). This way barriers of the order of kT or smaller do not hinder the random walk and a system started far far away from the important region in configuration space (blue lake) can move there quickly.

Importance sampling can be generalized to trajectory space, as we have done in our group. Consider the ensemble of all trajectories that are, say, 1 ps long. Most of these trajectories will be localized near some basin of attraction -- a long lived collection of neighboring micro states. Rare transition state crossings will comprise a small subset of these 1 ps trajectories. For example, if the process of interest occurs roughly once every milli second, then only one out of a billion 1 ps trajectories will exemplify that process. The method we have devised, called transition path sampling, harvests members of such rare sub ensembles.

 

Fig. 3: In a Transition Path Sampling simulation the path ensemble is sampled by generating a random walk in path space. This figure shows two steps of a such a random walk. Here, the blue path is generated by changing the red path by a small amount. Since the blue path is still connecting A with B it is accepted as the new path. By the same procedure the green path is generated from the blue path. This time, however, the new path is not a reactive trajectory. Therefore it is rejected. Indeed, it is this sequence of acceptances and rejections which ensures that the correct path ensemble is sampled in our computer simulation.

Suppose the rare processes of interest are transitions between states A and B, where A and B are characterized by their respective population operators, h_A(x) and h_B(x), as shown in Fig. 3. Here, x denotes a point in phase space. The probability weight for such unusual trajectories is

where p(x₀) is the unconstrained distribution of initial phase space points, x₀, and x_t is the phase space point t time steps later. (For simplicity, but not for necessity, we might be considering deterministic dynamics, in which case x_t is determined by the initial phase space point, x₀.) In this definition the population operators h_A(x) and h_B(x) are used to pick the reactive trajectories from the set of all possible trajectories. Therefore, the distribution P(x) is a statistical description of all trajectories going from A to B. We call this set of reactive paths the Transition Path Ensemble. Transition path sampling is done by carrying out a random walk in trajectory space, biased to be the importance sampling for the path distribution P(x).

In this perspective, one needs only to characterize stable or long lived states. For instance, one must characterize, perhaps by listing ranges of intramolecular lengths, an acetic acid molecule in water vs acetate and hydronium ions in water. Nothing is declared about the dynamical pathways (i.e., trajectories) that join these states. Rather, the transition path sampling is a random walk through the sub ensemble of all such possible paths. With this perspective, the calculation of time correlation functions becomes isomorphic with the calculation of reversible work or free energy differences [6]. Thus, standard tools of equilibrium computer simulation, like thermodynamic perturbation theory, can now be adopted to compute dynamical properties, like rate constants.

During the last year, we have demonstrated this approach in a series of published and soon to be published illustrations. These include studies of isomerizations in cold or low energy clusters of argon, and of water, and of ionic dissociation in water. The method can be melded with any procedure that generates reversible trajectories. For example, the method has been combined Car and Parrinello's ab initio molecular dynamics procedure in order to study the dissociation of weak acids in water. Transition path sampling has also been combined with Karplus' CHARMM program for atomistic modeling bio polymers.

In what follows on these pages, several published applications [4,5,7,8] of transition path sampling are briefly discussed. Reviews, one general [10], and one detailed [11] have been written. An up-to-date web page on transition path sampling is found here.

Applications

Even simple clusters of a few atoms can show interesting dynamical behavior. We studied rearrangement processes in a cluster consisting of seven Lennard-Jones particles in two dimensions. Using a quenching technique in path space we were able to identify the transition mechanisms without prior knowledge of the transition states. An isomorphism between time correlation functions and free energies can be employed to calculate rate constants.

Fig. 4: Potential energy profile obtained by transition path sampling for the isomerization of a 7-particle Lennard-Jones cluster. In the course of the transition the red particle migrates to its surface and is replaced by the blue particle at the center of the cluster. The configurations depicted below the curves are the stable states visited during the transition and the configurations above the curves are the transition states.

We have investigated proton transfer dynamics of the a protonated water trimer using transition path sampling. Transfer of a proton from one water molecule to another involves significant rearrangement of the cluster. Fig. 5 shows snapshots along a typical transition path.

 

Fig. 5: Snapshots of the protonated water trimer along a transition trajectory. During the transition a proton is transferred from water molecule 1 to water molecule 2 as molecule 3 moves from the left to the right side of the cluster.

During the transition, water molecule 3, which initially accepts a hydrogen bond from molecule 1, moves to the opposite side of the cluster passing over a saddle point on the potential energy surface. This rearrangement of the water molecules is accompanied by a proton transfer from molecule 1 to molecule 2. Two different transition state regions are important for the transition dynamics. They are shown in Fig. 6.

 

Fig. 6: View of the two relevant saddle points found using the path sampling procedure. In the higher energy transition state (left hand side)) the out-of-plane hydrogen atoms of water molecules 1 and 2 lie on the same side of the oxygen plane. In the lower energy transition state (right hand side) these hydrogens lie on opposite sides of the plane, so that the dipoles of water molecules 1 and 2 are antiparallel.

Since the transferring proton is a light particle, one might suspect that quantum mechanical effects are important for the transfer of the proton from one water to another. More specifically, one could imagine that to move from the donating water molecule to the accepting water molecule the proton has to cross a potential energy barrier. In this case, according to the bizarre rules of quantum mechanics, the proton can ''tunnel'' through the barrier momentarily violating the law of conservation of energy.

By using the path sampling procedure we discovered that the situation is different in the protonated water cluster. As mentioned above, proton transfer in this system is driven by the rearrangement of the water molecules. As a consequence the transferring proton feels a monostable potential trough the whole transition. Instead of tunneling through the barrier the proton is transported by the moving potential from one water to the other.

Imagine you put a pair of ions , say a Na⁺ ion and a Cl^- ion, in a glass of water. Due to the opposite charges they carry the two ions feel attracted to each other and tend to stay together. Eventually, however, thermal motion will lead to a separation of the two ions which then diffuse around freely in the liquid independently from each other. It turns out that in order to separate the two ions have to surmount a free energy barrier. Now, what is the role of the water molecules surrounding the ions in this barrier crossing process? Do the water molecules just provide the thermal energy to escape over the barrier or do they have a more specific role related to their molecular nature?

Fig. 8: A Na+Cl- ion pair immersed in water. Only the water molecules close to the ion pair are shown. Such an ion pair will dissociate with a characteristic time of about 15 picoseconds (1 picosecond = 10-12 seconds).

In room temperature water, a given water molecule dissociates on average once every 11 hours (a time scale nearly 20 orders of magnitude longer than the time scale of atomic motions!). Although the equilibrium concentration of the resulting species, aqueous hydronium (H₃O⁺) and hydroxide (OH^-) ions, is extremely low, the autodissociation process is fundamental to the chemistry of acids and bases in aqueous solution.

 

Fig. 9: Transition state for water autodissociation. The two illustrated configurations, separated by 30 femtoseconds(1 femtosecond = 10-15 seconds), show 32 water molecules shortly before (left) and after (right) crossing the transition state surface. The 32 molecules form the primary cell that is periodically replicated in space. The system is at room temperature, and obeys Newtonian dynamics with a force field determined by the BLYP electronic density functional theory, as determined by the Car-Parrinello molecular dynamics algorithm. The oxygen atoms colored yellow and blue highlight the hydronium and hydroxide ions, respectively. The dotted lines show hydrogen bonding proton wires along which ionic species move relatively quickly.

In the absence of a nearby counter ion, the Grotthus diffusion process occurs on a picosecond time scale. On the other hand, the presence of a nearby counter ion produces a strongly biased electric field that acts to force the diffusion process on a much more rapid time scale, and most typically in the direction of recombination. Temporal fluctuations of the solvent, however, can oppose recombination and indeed lead to the initial separation of hydroxide and hydronium ions. Transition states for autodissociation therefore involve the coincidence of two events. First, there must be a solvent fluctuation in local potential energy propelling the proton away from a nascent hydroxide ion. Second, this fluctuation must be accompanied by a breaking of a hydrogen bond along the proton wire connecting the hydroxide and hydronium ions. Without this second event, recombination will occur on the typically sub-picosecond time scale of the solvent fluctuation.

References

[1] This two step procedure was pioneered by Charles Bennett and David Chandler. See, for example, C. H. Bennett, in Algorithms for Chemical Computations, ACS Symp. Ser. No. 46, ed. R. E. Christofferson, American Chemical Society, Washington, 1977, p.63, and D. Chandler, J. Chem. Phys. 68, 2959 (1978). Related and earlier work has been discussed in a brief history communicated by D. Chandler, Faraday Discuss. Chem. Soc. 85, 341 (1988).

[2] Chandler, D., "Finding Transition Pathways: Throwing Ropes Over Rough Mountain Passes, in the Dark," in Computer Simulation of Rare Events and Dynamics of Classical and Quantum Condensed-Phase Systems -- Classical and Quantum Dynamics in Condensed Phase Simulations , edited by B. J. Berne, G. Ciccotti and D. F. Coker (World Scientific, Singapore, 1998), p 51-66. [pdf]

[3] C. Dellago, P. Bolhuis, F. S. Csajka, and D. Chandler, "Transition Path Sampling and the Calculation of Rate Constants", J. Chem. Phys., 108, 1964 (1998).[PDF]

[4] C. Dellago, P. Bolhuis, and D. Chandler, "Efficient Transition Path Sampling: Applications to Lennard-Jones cluster rearrangements", J. Chem. Phys., 108, 9236 (1998).[PDF]

[5] P. Bolhuis, C. Dellago and D. Chandler, "Sampling ensembles of deterministic transition pathways", Faraday Discuss., 110, 421 (1998).

[6] C. Dellago, P. G. Bolhuis and D. Chandler, "On the calculation of rate constants in the transition path ensemble, J. Chem. Phys., 110, 6617 (1999). [PDF]

[7] P. Geissler, C. Dellago, and D. Chandler, "Chemical dynamics of the protonated water trimer", Phys. Chem. Chem. Phys., 1, 1317 (1999). [PDF]

[8] P. Geissler, C. Dellago, and D. Chandler, "Kinetic Pathways of Ion Pair Dissociation in Water", J. Phys. Chem. B, 103, 3706 (1999). [PDF]

[9] Geissler, P.L., C. Dellago, D. Chandler, J. Hutter and M. Parinello, "Autoionization in liquid water", Science, 291, 212 (2001).[PDF]

[10]Bolhuis, P. G., D. Chandler, C. Dellago, and P. Geissler, "Transition Path Sampling: Throwing ropes over mountain passes, in the dark," Ann. Rev. of Phys. Chem., 59, 291-318 (2002). [PDF]

[11]C. Dellago, P. G. Bolhuis and P. L. Geissler, "Transition Path Sampling", Advances in Chemical Physics, 123, 1-78 (2002).[PDF]