Discrete path sampling

Multiscale dynamics of cluster fragmentation

The fragmentation of rare-gas clusters is theoretically investigated over time scales ranging from ionization and electronic excitation ( approximately fs) up to experimentally relevant times ( approximately ms). For this purpose a combination of methods are used, including nonadiabatic molecular dynamics, classical dynamics on the ground electronic state surface, and a kinetic description for the final evaporative cascade. The present multiscale protocol shows that, although the clusters are strongly out of equilibrium upon excitation, the long-time properties appear as statistical already after 1 ps.

Graph transformation method for calculating waiting times in Markov chains

We describe an exact approach for calculating transition probabilities and waiting times in finite-state discrete-time Markov processes. All the states and the rules for transitions between them must be known in advance. We can then calculate averages over a given ensemble of paths for both additive and multiplicative properties in a nonstochastic and noniterative fashion. In particular, we can calculate the mean first-passage time between arbitrary groups of stationary points for discrete path sampling databases, and hence extract phenomenological rate constants. We present a number of examples to demonstrate the efficiency and robustness of this approach.

Global optimization and folding pathways of selected alpha-helical proteins

The results of basin-hopping global optimization simulations are presented for four small, alpha-helical proteins described by a coarse-grained potential. A step-taking scheme that incorporates the local conformational preferences extracted from a large number of high-resolution protein structures is compared with an unbiased scheme. In addition, the discrete path sampling method is used to investigate the folding of one of the proteins, namely, the villin headpiece subdomain. Folding times from kinetic Monte Carlo simulations and iterative calculations based on a Markovian first-step analysis for the resulting stationary-point database are in good mutual agreement, but differ significantly from the experimental values, probably because the native state is not the global free energy minimum for the potential employed.

Sampling diffusive transition paths

The authors address the problem of sampling double-ended diffusive paths. The ensemble of paths is expressed using a symmetric version of the Onsager-Machlup formula, which only requires evaluation of the force field and which, upon direct time discretization, gives rise to a symmetric integrator that is accurate to second order. Efficiently sampling this ensemble requires avoiding the well-known stiffness problem associated with the sampling of infinitesimal Brownian increments of the path, as well as a different type of stiffness associated with the sampling of the coarse features of long paths. The fine-feature sampling stiffness is eliminated with the use of the fast sampling algorithm, and the coarse-feature sampling stiffness is avoided by introducing the sliding and sampling (S&S) algorithm. A key feature of the S&S algorithm is that it enables massively parallel computers to sample diffusive trajectories that are long in time. The authors use the algorithm to sample the transition path ensemble for the structural interconversion of the 38-atom Lennard-Jones cluster at low temperature.

Energy landscapes for water clusters in a uniform electric field

The behavior of small water clusters, (H2O)n, n=2-5 and n=8, in a uniform electric field is investigated for three related rigid-body models. Changes in the properties of the low-lying potential energy minima and the rearrangement pathways between them are examined. Results for certain structural transitions are compared with recent ab initio calculations. The models are found to give qualitatively similar trends, and there is some evidence that as the applied field strength is increased the quantitative differences between the models are also reduced.

Isomerization dynamics and thermodynamics of ionic argon clusters

The dynamics and thermodynamics of small Ar(n) (+) clusters, n=3, 6, and 9, are investigated using molecular dynamics (MD) and exchange Monte Carlo (MC) simulations. A diatomic-in-molecule Hamiltonian provides an accurate model for the electronic ground state potential energy surface. The microcanonical caloric curves calculated from MD and MC methods are shown to agree with each other, provided that the rigorous conservation of angular momentum is accounted for in the phase space density of the MC simulations. The previously proposed projective partition of the kinetic energy is used to assist MD simulations in interpreting the cluster dynamics in terms of inertial, internal, and external modes. The thermal behavior is correlated with the nature of the charged core in the cluster by computing a dedicated charge localization order parameter. We also perform systematic quenches to establish a connection with the various isomers. We find that the Ar(3) (+) cluster is very stable in its linear ground state geometry up to about 300 K, and then isomerizes to a T-shaped isomer in which a quasineutral atom lies around a charged dimer. In Ar(6) (+) and Ar(9) (+), the covalent trimer core is solvated by neutral atoms, and the weakly bound solvent shell melts at much lower energies, occasionally leading to a tetramer or pentamer core with weakly charged extremities. At high energies the core itself becomes metastable and the cluster transforms into Ar(2) (+) solvated by a fluid of neutral argon atoms.

Relaxation dynamics and structural isomerism in Nb10 and Nb10+

The structure, energetics, and interconversion of isomers of Nb10 and Nb10+ are studied using density functional theory with Gaussian basis sets, using guess structures derived from basin-hopping simulations with the Finnis-Sinclair [Philos. Mag. A 50, 45 (1984)] potential. These results are used as input to a master equation approach to model the relaxation of these clusters. Ionization potentials are calculated for all relevant minima, as are the infrared spectra. On the basis of these data, and known experimental results, plausible explanations are given for the biexponential reaction kinetics observed for Nb10 and Nb10+ with respect to small molecule adsorbates. In principle, this approach could be extended to investigate any midsized transition metal cluster that exhibits structural isomerism.

A doubly nudged elastic band method for finding transition states

A modification of the nudged elastic band (NEB) method is presented that enables stable optimizations to be run using both the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) quasi-Newton and slow-response quenched velocity Verlet minimizers. The performance of this new "doubly nudged" DNEB method is analyzed in conjunction with both minimizers and compared with previous NEB formulations. We find that the fastest DNEB approach (DNEB/L-BFGS) can be quicker by up to 2 orders of magnitude. Applications to permutational rearrangements of the seven-atom Lennard-Jones cluster (LJ7) and highly cooperative rearrangements of LJ38 and LJ75 are presented. We also outline an updated algorithm for constructing complicated multi-step pathways using successive DNEB runs.

Transition rate prefactors for systems of many degrees of freedom

When a minimum on the potential energy surface is surrounded by multiple saddle points with similar energy barriers, the transition pathways with greater prefactors are more important than those that have similar energy barriers but smaller prefactors. In this paper, we present a theoretical formulation for the prefactors, computing the probabilities for transition paths from a minimum to its surrounding saddle points. We apply this formulation to a system of 2 degrees of freedom and a system of 14 degrees of freedom. The first is Brownian motion in a two-dimensional potential whose global anharmonicities play a dominant role in determining the transition rates. The second is a Lennard-Jones (LJ) cluster of seven particles in two dimensions. Low lying transition states of the LJ cluster, which can be reached directly from a minimum without passing through another minimum, are identified without any presumption of their characteristics nor of the product states they lead to. The probabilities are computed for paths going from an equilibrium ensemble of states near a given minimum to the surrounding transition states. These probabilities are directly related to the prefactors in the rate formula. This determination of the rate prefactors includes all anharmonicities, near or far from transition states, which are pertinent in the very sophisticated energy landscape of LJ clusters and in many other complex systems.

Energy landscapes and properties of biomolecules

Thermodynamic and dynamic properties of biomolecules can be calculated using a coarse-grained approach based upon sampling stationary points of the underlying potential energy surface. The superposition approximation provides an overall partition function as a sum of contributions from the local minima, and hence functions such as internal energy, entropy, free energy and the heat capacity. To obtain rates we must also sample transition states that link the local minima, and the discrete path sampling method provides a systematic means to achieve this goal. A coarse-grained picture is also helpful in locating the global minimum using the basin-hopping approach. Here we can exploit a fictitious dynamics between the basins of attraction of local minima, since the objective is to find the lowest minimum, rather than to reproduce the thermodynamics or dynamics.

Finding pathways between distant local minima

We report a new algorithm for constructing pathways between local minima that involve a large number of intervening transition states on the potential energy surface. A significant improvement in efficiency has been achieved by changing the strategy for choosing successive pairs of local minima that serve as endpoints for the next search. We employ Dijkstra's algorithm [E. W. Dijkstra, Numer. Math. 1, 269 (1959)] to identify the "shortest" path corresponding to missing connections within an evolving database of local minima and the transition states that connect them. The metric employed to determine the shortest missing connection is a function of the minimized Euclidean distance. We present applications to the formation of buckminsterfullerene and to the folding of various biomolecules: the B1 domain of protein G, tryptophan zippers, and the villin headpiece subdomain. The corresponding pathways contain up to 163 transition states and will be used in future discrete path sampling calculations.

On the Convergence of Parallel Tempering Monte Carlo Simulations of LJ38

The convergence of parallel tempering Monte Carlo simulations of the 38-atom Lennard-Jones cluster starting from the Oh global minimum and from the C(5v) second-lowest-energy minimum is investigated. It is found that achieving convergence is appreciably more difficult, particularly at temperatures in the vicinity of the Oh --> C(5v) transformation when starting from the C(5v) structure. A strategy combining the Tsallis generalized ensemble and the parallel tempering algorithm is implemented and used to improve the convergence of the simulations in the vicinity of the Oh --> C(5v) transformation.

Protonated water clusters described by an empirical valence bond potential

The properties of low-lying stationary points on the potential energy surfaces of singly protonated water clusters (H(2)O)(n)H(+), are investigated using an empirical valence bond potential. Candidate global minima are reported for n=2-4, 8, and 20-22. For n=8, the variation in the energies and structures of low-lying minima with the number of valence bond states included in the model is studied. For n=4 and 8, disconnectivity graphs are also reported and are compared to results for the equivalent neutral water clusters as described by the rigid TIP3P potential. For the larger clusters, n=20-22, the structural properties of the low energy minima are compared with recently published spectroscopic data on these systems. The observed differences between the n=20 and n=21 systems are qualitatively reproduced by the model potential, but the similarities between the n=21 and n=22 systems are not.

Free energy profile along a discretized reaction path via the hyperplane constraint force and torque

By employing mechanical work analogies, we derive a convenient computational approach for evaluation of the free energy profile (FEP) along some discretized path defined as a sequence of hyperplanes. A hyperplane is fully specified by any of its point and a tangent vector. The FEP is obtained as an integral of two components. The translational component of the free energy is computed by integrating the hyperplane constraint force. The rotational component is evaluated via the hyperplane torque. Both ingredients--the constraint force and the hyperplane torque-are evaluated on each hyperplane independently. The integration procedure utilizes a set of reference points defining a point of rotation on each hyperplane, and these points can be chosen before or after the sampling takes place. A shift in the reference points redistributes the FEP contributions between the translational and rotational components. For systems where the FEP is dominated by the potential energy differences, reference points residing on the minimum energy path present a natural choice. We demonstrate the validity of our approach on two examples, a simple two-dimensional (2D) potential, and a seven-atom Lennard-Jones cluster. In each case, we compare the numerical FEP with the harmonic approximation estimates. Our results for the 2D potential are also verified by the data available in the literature. In both cases, the rotational component of the FEP represents a sizable contribution to the total FEP, so ignoring it would yield clearly incorrect results.

Characterizing the network topology of the energy landscapes of atomic clusters

By dividing potential energy landscapes into basins of attractions surrounding minima and linking those basins that are connected by transition state valleys, a network description of energy landscapes naturally arises. These networks are characterized in detail for a series of small Lennard-Jones clusters and show behavior characteristic of small-world and scale-free networks. However, unlike many such networks, this topology cannot reflect the rules governing the dynamics of network growth, because they are static spatial networks. Instead, the heterogeneity in the networks stems from differences in the potential energy of the minima, and hence the hyperareas of their associated basins of attraction. The low-energy minima with large basins of attraction act as hubs in the network. Comparisons to randomized networks with the same degree distribution reveals structuring in the networks that reflects their spatial embedding.

The energy landscape as a unifying theme in molecular science

The potential energy surface (PES) underlies most calculations of structure, dynamics and thermodynamics in molecular science. In this contribution connections between the topology of the PES and observable properties are developed for a coarse-grained model of virus capsid self-assembly. The model predicts that a thermodynamically stable, kinetically accessible icosahedral shell exists for pentameric building blocks of the right shape: not too flat and not too spiky. The structure of the corresponding PES is probably common to other systems where directed searches avoid Levinthal's paradox, such as 'magic number' clusters, protein folding and crystallization.

Archetypal energy landscapes: Dynamical diagnosis

Recent studies have identified several motifs for potential energy surfaces corresponding to distinct dynamic and thermodynamic properties. The corresponding disconnectivity graphs were identified as "palm tree," "willow tree," and "banyan tree" patterns. In the present contribution we present a quantitative analysis of the relation between the topography and dynamics for each of these motifs. For the palm tree and willow tree forms we find that the arrangement of the stationary points in the monotonic sequences with respect to the global minimum is the most important factor in establishing the kinetic properties. However, the results are somewhat different for motifs involving a rough surface with several deep basins (banyan tree motif), with large barriers relative to the energy differences between minima. Here it is the size of the barrier for escape from the region relative to the barriers at the bottom that is most important. The present results may be helpful in distinguishing between the dynamics of "structure seeking" and "glass forming" systems.

Tools for channels: moving towards molecular calculations of gating and permeation in ion channel biophysics

Recent X-ray structures of voltage gated potassium channels provide an exciting opportunity to connect molecular structures with measured biological function. Two of the most important connections for these channels are: first, to the molecular basis behind selectivity and the associated free energy profile underlying ionic current flow and, second, to a true molecular understanding of the large-scale conformational transitions that underlie voltage dependent gating. But, existing computational tools need to be further developed to reach these goals. In this contribution to the symposia on sampling methods we outline our dynamic importance sampling method for sampling large-scale conformational transitions as well as our studies with non-equilibrium work events and equilibrium overlap sampling (OS) methods for sampling events related to the calculation of relative free energies.

Analysis of cooperativity and localization for atomic rearrangements

We propose measures of localization and cooperativity for the analysis of atomic rearrangements. We show that for both clusters and bulk material cooperative rearrangements usually have significantly lower barriers than uncooperative ones, irrespective of the degree of localization. We also find that previous methods used to sample stationary points are biased towards rearrangements of particular types. Linear interpolation between local minima in double-ended transition state searches tends to produce cooperative rearrangements, while random perturbations of all the coordinates, as sometimes used in single-ended searches, have the opposite effect.

Hidden complexity of free energy surfaces for peptide (protein) folding

An understanding of the thermodynamics and kinetics of protein folding requires a knowledge of the free energy surface governing the motion of the polypeptide chain. Because of the many degrees of freedom involved, surfaces projected on only one or two progress variables are generally used in descriptions of the folding reaction. Such projections result in relatively smooth surfaces, but they could mask the complexity of the unprojected surface. Here we introduce an approach to determine the actual (unprojected) free energy surface and apply it to the second beta-hairpin of protein G, which has been used as a model system for protein folding. The surface is represented by a disconnectivity graph calculated from a long equilibrium folding-unfolding trajectory. The denatured state is found to have multiple low free energy basins. Nevertheless, the peptide shows exponential kinetics in folding to the native basin. Projected surfaces obtained from the present analysis have a simple form in agreement with other studies of the beta-hairpin. The hidden complexity found for the beta-hairpin surface suggests that the standard funnel picture of protein folding should be revisited.

Folding of the GB1 hairpin peptide from discrete path sampling

The discrete path sampling technique is used to calculate folding pathways of the 16-amino acid beta hairpin-forming sequence from residues 41-56 of the B1 domain of protein G. The folding time is obtained using master equation dynamics and kinetic Monte Carlo simulations, and the time evolution of different order parameters and occupation probabilities of groups of minima are calculated and used to characterize intermediates on the folding pathway.

Comparison of kinetic Monte Carlo and molecular dynamics simulations of diffusion in a model glass former

We present results from kinetic Monte Carlo (KMC) simulations of diffusion in a model glass former. We find that the diffusion constants obtained from KMC simulations have Arrhenius temperature dependence, while the correct behavior, obtained from molecular dynamics simulations, can be super-Arrhenius. We conclude that the discrepancy is due to undersampling of higher-lying local minima in the KMC runs. We suggest that the relevant connectivity of minima on the potential energy surface is proportional to the energy density of the local minima, which determines the "inherent structure entropy." The changing connectivity with potential energy may produce a correlation between dynamics and thermodynamics.