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.     

Steered transition path sampling

We introduce a path sampling method for obtaining statistical properties of an arbitrary stochastic dynamics. The method works by decomposing a trajectory in time, estimating the probability of satisfying a progress constraint, modifying the dynamics based on that probability, and then reweighting to calculate averages. Because the progress constraint can be formulated in terms of occurrences of events within time intervals, the method is particularly well suited for controlling the sampling of currents of dynamic events. We demonstrate the method for calculating transition probabilities in barrier crossing problems and survival probabilities in strongly diffusive systems with absorbing states, which are difficult to treat by shooting. We discuss the relation of the algorithm to other methods.     

Comparison and union of the Temple and Bazley lower bounds

The Lehmann-Maehly approach and Bazley’s method of special choice are matrix eigenvalue problems that allow the calculation of lower bounds to energies of atomic and molecular systems. We introduce a common derivation of their scalar versions using the overlap of a trial function with the unknown ground-state wave function. In the scalar setting, the Lehmann-Maehly approach reduces to the Temple formula. The common derivation allows us to easily unite and improve both methods in several stages within this restricted application. Finally we offer a different union that allows generalization to arbitrary dimension matrix methods. Calculations on the helium atom ground state illustrate the improvements and mergers.     

The role of model systems in the few-body reduction of the N-fermion problem

The problem of upper and lower ground state energy bounds for many-fermion systems is considered from the viewpoint of reduced density matrices. Model density matrices are used for upper bounds to, first uncoupled, then coupled fermions. Model Hamiltonians are developed for lower bounds in corresponding fashion. Both mathematical and physical models are constructed for setting up universally valid inequalities on density matrices. These are joined by both inequalities and equalities in which the explicit form of the system at hand is used. A few illustrative examples are presented.     

Theory of photoinduced heterogeneous electron transfer

We consider electron injection into the conduction band of a semiconductor, from an electronically excited state of a dye molecule, adsorbed on its surface. For arbitrary width of the conduction band, the survival probability of the excited state can be calculated using a Green's-function approach. We show that the existence of a split-off state can play an important role in the total injection probability. In the wide band limit, the survival probability decays exponentially, but for finite band widths it does not. We further investigate the effect of vibrations on the process. A Green's operator technique may be used to solve this too exactly. We show that the problem may be reduced to a non-Hermitian eigenvalue problem for the vibrational states alone. Exact results can be obtained for arbitrary bandwidth and for a few vibrational degrees of freedom. In the wide band limit, the dynamics is particularly simple and we find that (1) the survival probability of the excited state is unchanged by the inclusion of vibrational motion, but (2) each vibrational state now has a finite lifetime. Numerical results are presented for the effects of reorganization energy, energy of the injecting level, and the variation of the matrix element for the electron injection, on the survival probability of the electron in the excited state. As an illustration of the approach, we also present results of numerical calculation of the absorption spectrum of perylene adsorbed on TiO(2) and compare it with experimental results.     

Dependence of the recombination kinetics of geminate electron-ion pairs on the mean free path of electrons

We have investigated how the kinetics of geminate recombination of electron-ion pairs is influenced by the mean free path of electrons. By use of the Monte Carlo method previously developed the decay of the pair survival probability was numerically calculated for a variety of values of the mean free path of electrons. When the mean free path is negligibly small compared with the Onsager length, our result is in good agreement with the Hong-Noolandi result based on the Smoluchowski equation. However, as the mean free path increases, it deviates from their result increasingly. This indicates that the treatment based on the Smoluchowski equation is inadequate when the mean free path of electrons is not negligible compared with the Onsager length. We have found that for mean free paths comparable with or larger than the Onsager length the decay kinetics of the pair survival probability approximately follows the first-order kinetics except at very short and very long times. The dependence of the escape probability on the mean free path of electrons is also briefly discussed.     

Accurate free energy calculation along optimized paths

The path‐based methods of free energy calculation, such as thermodynamic integration and free energy perturbation, are simple in theory, but difficult in practice because in most cases smooth paths do not exist, especially for large molecules. In this article, we present a novel method to build the transition path of a peptide. We use harmonic potentials to restrain its nonhydrogen atom dihedrals in the initial state and set the equilibrium angles of the potentials as those in the final state. Through a series of steps of geometrical optimization, we can construct a smooth and short path from the initial state to the final state. This path can be used to calculate free energy difference. To validate this method, we apply it to a small 10‐ALA peptide and find that the calculated free energy changes in helix‐helix and helix‐hairpin transitions are both self‐convergent and cross‐convergent. We also calculate the free energy differences between different stable states of β‐hairpin trpzip2, and the results show that this method is more efficient than the conventional molecular dynamics method in accurate free energy calculation. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Adaptive importance sampling Monte Carlo simulation of rare transition events

We develop a general theoretical framework for the recently proposed importance sampling method for enhancing the efficiency of rare-event simulations [W. Cai, M. H. Kalos, M. de Koning, and V. V. Bulatov, Phys. Rev. E 66, 046703 (2002)], and discuss practical aspects of its application. We define the success/fail ensemble of all possible successful and failed transition paths of any duration and demonstrate that in this formulation the rare-event problem can be interpreted as a "hit-or-miss" Monte Carlo quadrature calculation of a path integral. The fact that the integrand contributes significantly only for a very tiny fraction of all possible paths then naturally leads to a "standard" importance sampling approach to Monte Carlo (MC) quadrature and the existence of an optimal importance function. In addition to showing that the approach is general and expected to be applicable beyond the realm of Markovian path simulations, for which the method was originally proposed, the formulation reveals a conceptual analogy with the variational MC (VMC) method. The search for the optimal importance function in the former is analogous to finding the ground-state wave function in the latter. In two model problems we discuss practical aspects of finding a suitable approximation for the optimal importance function. For this purpose we follow the strategy that is typically adopted in VMC calculations: the selection of a trial functional form for the optimal importance function, followed by the optimization of its adjustable parameters. The latter is accomplished by means of an adaptive optimization procedure based on a combination of steepest-descent and genetic algorithms.     

Higher order Morita approximations for random copolymer localization

The Morita approximation is a constrained annealing procedure which yields upper bounds on the quenched average free energy for models of quenched randomness. In this article we consider a bilateral Dyck path model, first introduced by S.G. Whittington and collaborators, of the localization of a random copolymer at the interface between two immiscible solvents. The distribution of comonomers along the polymer chain is initially determined by a random process and once chosen it remains fixed. Morita approximations in which we control correlations to various orders between neighbouring monomers along the polymer chain are applied to this model. Although at low orders the Morita approximation does not yield the correct path properties in the localized region of the phase diagram, we show that this problem can be overcome by including sufficiently high-order correlations in the Morita approximation. In addition by comparison with an appropriate lower bound, we show that well-within the localized phase the Morita approximation provides a relatively tight upper bound on the limiting quenched average free energy for bilateral Dyck path localization.     

Activation processes with memory

We propose a mathematical treatment of the activated processes governed by stochastic Langevin dynamics with a colored random force, corresponding to a noise generated by an Ornstein-Uhlenbeck process. Such non-Markovian dynamics take place in a variety of chemical and biological systems. Using the path integral approach, we constructed the conditional probability for passing between two stationary states in configurational space. Our relations can be used for Monte Carlo sampling of evolution trajectories for systems with many degrees of freedom as well as for determining the reaction coordinate used in transition state theory. On the basis of our relation for a conditional probability, we generalize the method of determining the most probable path to the case of colored random force. Using the simple three-hole potential, we examine numerically the effect of nonzero correlation time (memory) on the evolution of the most probable path for a finite temperature.     

Modeling flexible pharmacophores with distance geometry, scoring, and bound stretching

The study of pharmacophores, i.e., of common features between different ligands, is important for the quantitative identification of "compatible" enzymes and binding species. A pharmacophore-based technique is developed that combines multiple conformations with a distance geometry method to create flexible pharmacophore representations. It uses a set of low-energy conformations combined with a new process we call bound stretching to create sets of distance bounds, which contain all or most of the low-energy conformations. The bounds can be obtained using the exact distances between pairs of atoms from the different low-energy conformations. To avoid missing conformations, we can take advantage of the triangle distance inequality between sets of three points to logically expand a set of upper and lower distance bounds (bound stretching). The flexible pharmacophore can be found using a 3-D maximal common subgraph method, which uses the overlap of distance bounds to determine the overlapping structure. A scoring routine is implemented to select the substructures with the largest overlap because there will typically be many overlaps with the maximum number of overlapping bounds. A case study is presented in which 3-D flexible pharmacophores are generated and used to eliminate potential binding species identified by a 2-D pharmacophore method. A second case study creates flexible pharmacophores from a set of thrombin ligands. These are used to compare the new method with existing pharmacophore identification software.     

Use of augmented Lagrangians in the calculation of molecular conformations by distance geometry

Distance geometry is a technique widely used to find atomic coordinates that agree with given upper and lower bounds on the interatomic distances. It is successful because it chooses at random some relatively good "trial coordinates" that take into account the whole molecule and all constraints at once. Customarily, these trial coordinates must be refined by minimizing a penalty function until the structure agrees with the original bounds. Here we present an alternative to minimizing the penalty function, which has the advantage of more precisely satisfying the bounds, showing more clearly when the bounds are mutually contradictory, and simultaneously optimizing an objective function subject to precise satisfaction of the bounds.     

Maximum Flux Transition Paths of Conformational Change

Given two metastable states A and B of a biomolecular system, the problem is to calculate the likely paths of the transition from A to B. Such a calculation is more informative and more manageable if done for a reduced set of collective variables chosen so that paths cluster in collective variable space. The computational task becomes that of computing the "center" of such a cluster. A good way to define the center employs the concept of a committor, whose value at a point in collective variable space is the probability that a trajectory at that point will reach B before A. The committor "foliates" the transition region into a set of isocommittors. The maximum flux transition path is defined as a path that crosses each isocommittor at a point which (locally) has the highest crossing rate of distinct reactive trajectories. This path is based on the same principle as the minimum resistance path of Berkowitz et al (1983), but it has two advantages: (i) the path is invariant with respect to a change of coordinates in collective variable space and (ii) the differential equations that define the path are simpler. It is argued that such a path is nearer to an ideal path than others that have been proposed with the possible exception of the finite-temperature string method path. To make the calculation tractable, three approximations are introduced, yielding a path that is the solution of a nonsingular two-point boundary-value problem. For such a problem, one can construct a simple and robust algorithm. One such algorithm and its performance is discussed.     

An algorithm to find minimum free-energy paths using umbrella integration

The calculation of free-energy barriers by umbrella sampling and many other methods is hampered by the necessity for an a priori choice of the reaction coordinate along which to sample. We avoid this problem by providing a method to search for saddle points on the free-energy surface in many coordinates. The necessary gradients and Hessians of the free energy are obtained by multidimensional umbrella integration. We construct the minimum free-energy path by following the gradient down to minima on the free-energy surface. The change of free energy along the path is obtained by integrating out all coordinates orthogonal to the path. While we expect the method to be applicable to large systems, we test it on the alanine dipeptide in vacuum. The minima, transition states, and free-energy barriers agree well with those obtained previously with other methods.     

Bounds for Hartree-Fock perturbation theory

Upper and lower bounds for the second-order energy in both coupled and uncoupled Hartree-Fock perturbation theories are derived. Using these bounds inequalities are derived for the error in the geometric approximation.     

Calculation of the Entropy of Lattice Polymer Models from Monte Carlo Trajectories

While lattice models are used extensively for macromolecules (synthetic polymers proteins, etc.), calculation of the absolute entropy, S, and the free energy, F, from a given Monte Carlo (MC) trajectory is not straightforward. Recently, we have developed the hypothetical scanning MC (HSMC) method for calculating S and F of fluids. Here we extend HSMC to self-avoiding walks on a square lattice and discuss its wide applicability to complex polymer lattice models. HSMC is independent of existing techniques and thus constitutes an independent research tool; it provides rigorous upper and lower bounds for F, which can be obtained from a very small sample and even from a single chain conformation.     

Generalised technique for calculation of plane director profiles in bounded nematic liquid crystals

We propose an improved conformal mapping technique for analytical calculation of two-dimensional profiles of the nematic liquid crystal (NLC) director in a cell with a simply connected cross-sectional region. We consider the case of the strong anchoring and the piecewise constant director pretilt on the piecewise smooth curve which bounds the region. Obtained expressions for the director profile explicitly depend on the conformal mapping which maps the region onto the upper half plane of a complex plane. An advantage of our method in comparison with the standard conformal mapping technique is that it does not require the knowledge of the inverse mapping and the calculation of the integral in the Poisson formula. Proposed technique allows to take into account topological defects in the bulk of the NLC on the symmetry axis of the region. As an example of how the method can be used, we find an analytical expression for the director profile in a horizontal cylindrical groove partly filled with the NLC. We consider the case where a disclination line parallel to the axis of the groove occurs in the bulk of the NLC. The equilibrium position of the disclination line is found.     

Variational bounds to the overlap

Starting from a closed expression for the overlap, variational upper and lower bounds to the overlap are derived by means of operator inequalities.