A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives

Transition pathways in complex systems: Application of the finite-temperature string method to the alanine dipeptide

The finite-temperature string method proposed by E, et al. [W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B 66, 052301 (2002)] is a very effective way of identifying transition mechanisms and transition rates between metastable states in systems with complex energy landscapes. In this paper, we discuss the theoretical background and algorithmic details of the finite-temperature string method, as well as the application to the study of isomerization reaction of the alanine dipeptide, both in vacuum and in explicit solvent. We demonstrate that the method allows us to identify directly the isocommittor surfaces, which are approximated by hyperplanes, in the region of configuration space where the most probable transition trajectories are concentrated. These results are verified subsequently by computing directly the committor distribution on the hyperplanes that define the transition state region.

DFT Characterization of Adsorbed NHx Species on Pt(100) and Pt(111) Surfaces

Periodic density functional theory (DFT) calculations using plane waves have been performed to systematically investigate the adsorption and relative stability of ammonia and its dehydrogenated species on Pt(111) and Pt(100) surfaces. Different adsorption geometries and positions have been studied, and in each case, the equilibrium configuration has been determined by relaxation of the system. The vibrational spectra of the various ammonia fragments have been computed, and band assignments have been compared in detail with available experimental data. The adsorption of NH3 (on top) and NH2 (bridge) is more favorable on Pt(100) than on Pt(111), while similar adsorption energies were computed for NH (hollow) and N (hollow) on both surfaces. The remarkably lower adsorption energy of NH2 over Pt(111) as compared with Pt(100) (the difference being approximately 0.7 eV) can be related to different geometric and electronic factors associated with this particular intermediate. Accordingly, the type of platinum surface determines the most stable NH(x) fragment: Pt(100) has more affinity for NH2 species, whereas NH species are preferred over Pt(111).

Selectivity Control for the Catalytic 1,3-Butadiene Hydrogenation on Pt(111) and Pd(111) Surfaces: Radical versus Closed-Shell Intermediates

The hydrogenation of 1,3-butadiene to different C4H8 species on both Pd(111) and Pt(111) surfaces has been studied by means of periodic slabs and DFT. We report the adsorption structures for the various mono- and dihydrogenated butadiene intermediates adsorbed on both metal surfaces. Radical species are more clearly stabilized on Pt than on Pd. The different pathways leading to these radicals have been investigated and compared to those producing 1-butene and 2-butene species. On palladium, the formation of butenes seems to be clearly favored, in agreement with the high selectivity to butenes observed experimentally. In contrast, the formation of dihydrogenated radical species seems to be competitive with that of butenes on platinum, which could explain its poorer selectivity to butenes and the formation of butane as a primary product.

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.

Finding transition states for crystalline solid-solid phase transformations

We present a method to identify transition states and minimum energy paths for martensitic solid-solid phase transformations, thereby allowing quantification of the activation energies of such transformations. Our approach is a generalization of a previous method for identifying transition states for chemical reactions, namely the climbing image-nudged elastic band algorithm, where here the global deformation of the crystalline lattice (volume and shape fluctuations) becomes the reaction coordinate instead of atomic motion. We also introduce an analogue to the Born-Oppenheimer approximation that allows a decoupling of nuclear motion and lattice deformation, where the nuclear positions along the path are determined variationally according to current deformation state. We then apply this technique to characterize the energetics of elemental lithium phase transformations as a function of applied pressure, where we see a validation of the Born-Oppenheimer-like approximation, small energy barriers (expected for martensitic transformations), and a pronounced pressure dependence of various properties characterizing the phase transitions.

Reaction path determination for quantum mechanical/molecular mechanical modeling of enzyme reactions by combining first order and second order “chain-of-replicas” methods

A two-step procedure for the determination of reaction paths in enzyme systems is presented. This procedure combines two chain-of-states methods: a quantum mechanical/molecular mechanical (QM/MM) implementation of the nudged elastic band (NEB) method and a second order parallel path optimizer method both recently developed in our laboratory. In the first step, a reaction path determination is performed with the NEB method, along with a restrained minimization procedure for the MM environment to obtain a first approximation to the reaction path. In the second step, the calculated path is refined with the parallel path optimizer method. By combining these two methods the reaction paths are determined accurately, and in addition, the number of path optimization iterations are significantly reduced. This procedure is tested by calculating both steps of the isomerization of 2-oxo-4-hexenedioate by 4-oxalocrotonate tautomerase, which have been previously determined by our group. The calculated paths agree with the previously reported results and we obtain a reduction of 45%-55% in the number of path optimization cycles.

A Recipe for the Computation of the Free Energy Barrier and the Lowest Free Energy Path of Concerted Reactions

The recently introduced hills method (Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562) is a powerful tool to compute the multidimensional free energy surface of intrinsically concerted reactions. We have extended this method by focusing our attention on localizing the lowest free energy path that connects the stable reactant and product states. This path represents the most probable reaction mechanism, similar to the zero temperature intrinsic reaction coordinate, but also includes finite temperature effects. The transformation of the multidimensional problem to a one-dimensional reaction coordinate allows for accurate convergence of the free energy profile along the lowest free energy path using standard free energy methods. Here we apply the hills method, our lowest free energy path search algorithm, and umbrella sampling to the prototype S(N)2 reaction. The hills method replaces the in many cases difficult problem of finding a good reaction coordinate with choosing relatively simple collective variables, such as the bond lengths of the broken and formed chemical bonds. The second part of the paper presents a guide to using the hills method, in which we test and fine-tune the method for optimal accuracy and efficiency using the umbrella sampling results as a reference.

Potential energy constrained molecular dynamics simulations

A method for carrying out molecular dynamics simulations in which the potential energy U of the molecular system is constrained at its initial value is developed and thoroughly tested. The constraint is not introduced within the framework of the Lagrange multipliers technique, rather it is fulfilled in a natural way by carrying out the simulations in terms of suitable sets of delocalized coordinates. Such coordinates are defined by an appropriate tuning of the Baker, Kessi, and Delley internal delocalized nonredundant coordinates technique [J. Chem. Phys. 105, 192 (1996)]. The proposed method requires multiple evaluations of energy and gradients in each step of the molecular dynamics simulation, so that constant U simulations suffer some overhead compared to ordinary simulations. But the particular formulation of the delocalized coordinates and of the equations of motion greatly simplifies all the various steps required by the Baker's technique, thus allowing for the efficient implementation of the method itself. The technique is reliable and allows for very high accuracy in the potential energy conservation during the whole simulation. Moreover, it proved to be free of drift troubles which can occur when standard constraint methods are straightforwardly implemented without the application of appropriate correcting techniques.

Mapping potential energy surfaces

A recently proposed dynamical method [A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. 99, 12562 (2002)] allows us to globally sample the free energy surface. This approach uses a coarse-grained non-Markovian dynamics to bias microscopic atomic trajectories. After a sufficiently long simulation time, the global free energy surface can be reconstructed from the non-Markovian dynamics. Here we apply this scheme to study the T=0 free energy surface, i.e., the potential energy surface in coarse-grained space. We show that the accuracy of the reconstructed potential energy surface can be dramatically improved by a simple postprocessing procedure with only minor computational overhead. We illustrate this approach by conducting conformational analysis on a small organic molecule, demonstrating its superiority over traditional unbiased approaches in sampling potential energy surfaces in coarse-grained space.

Atomistic study of dislocation loop emission from a crack tip

We report the first atomistic calculation of the saddle-point configuration and activation energy for the nucleation of a 3D dislocation loop from a stressed crack tip in single crystal Cu. The transition state is found using reaction pathway sampling schemes, the nudged elastic band, and dimer methods. For the (111)[110] crack, loaded typically at 75% of the athermal critical strain energy release rate for spontaneous dislocation nucleation, the calculated activation energy is 1.1 eV, significantly higher than the continuum estimate. Implications concerning homogeneous dislocation nucleation in the presence of a crack-tip stress field are discussed.

Influence of confinement on the vibrational density of states and the Boson peak in a polymer glass

We have performed a normal-mode analysis on a glass forming polymer system for bulk and free-standing film geometries prepared under identical conditions. It is found that for free-standing film glasses, the normal-mode spectrum exhibits significant differences from the bulk glass with the appearance of an additional low-frequency peak and a higher intensity at the Boson peak frequency. A detailed eigenvector analysis shows that the low-frequency peak corresponds to a shear-horizontal mode which is predicted by continuum theory. The peak at higher frequency (Boson peak) corresponds to motions that are correlated over a length scale of approximately twice the interaction site diameter. These observations shed some light on the microscopic dynamics of glass formers, and help explain decreasing fragility that arises with decreasing thickness in thin films.

A growing string method for determining transition states: Comparison to the nudged elastic band and string methods

Interpolation methods such as the nudged elastic band and string methods are widely used for calculating minimum energy pathways and transition states for chemical reactions. Both methods require an initial guess for the reaction pathway. A poorly chosen initial guess can cause slow convergence, convergence to an incorrect pathway, or even failed electronic structure force calculations along the guessed pathway. This paper presents a growing string method that can find minimum energy pathways and transition states without the requirement of an initial guess for the pathway. The growing string begins as two string fragments, one associated with the reactants and the other with the products. Each string fragment is grown separately until the fragments converge. Once the two fragments join, the full string moves toward the minimum energy pathway according to the algorithm for the string method. This paper compares the growing string method to the string method and to the nudged elastic band method using the alanine dipeptide rearrangement as an example. In this example, for which the linearly interpolated guess is far from the minimum energy pathway, the growing string method finds the saddle point with significantly fewer electronic structure force calculations than the string method or the nudged elastic band method.

Structure and mobility of defects formed from collision cascades in MgO

We study radiation-damage events in MgO on experimental time scales by augmenting molecular dynamics cascade simulations with temperature accelerated dynamics, molecular statics, and density functional theory. At 400 eV, vacancies and mono- and di-interstitials form, but often annihilate within milliseconds. At 2 and 5 keV, larger clusters can form and persist. While vacancies are immobile, interstitials aggregate into clusters (In) with surprising properties; e.g., an I4 is immobile, but an impinging I2 can create a metastable I6 that diffuses on the nanosecond time scale but is stable for years.

Competitive Paths for Methanol Decomposition on Pt(111)

Periodic, self-consistent, Density Functional Theory (PW91-GGA) calculations are used to study competitive paths for the decomposition of methanol on Pt(111). Pathways proceeding through initial C-H and C-O bond scission events in methanol are considered, and the results are compared to data for a pathway proceeding through an initial O-H scission event [Greeley et al. J. Am. Chem. Soc. 2002, 124, 7193]. The DFT results suggest that methanol decomposition via CH(2)OH and either formaldehyde or HCOH intermediates is an energetically feasible pathway; O-H scission to CH(3)O, followed by sequential dehydrogenation, may be another realistic route. Microkinetic modeling based on the first-principles results shows that, under realistic reaction conditions, C-H scission in methanol is the initial decomposition step with the highest net rate. The elementary steps of all reaction pathways (with the exception of C-O scission) follow a linear correlation between the transition state and final state energies. Simulated HREELS spectra of the intermediates show good agreement with available experimental data, and HREELS spectra of experimentally elusive reaction intermediates are predicted.

Efficient exploration of reactive potential energy surfaces using Car-Parrinello molecular dynamics

The possibility of observing chemical reactions in ab initio molecular dynamics runs is severely hindered by the short simulation time accessible. We propose a new method for accelerating the reaction process, based on the ideas of the extended Lagrangian and coarse-grained non-Markovian metadynamics. We demonstrate that by this method it is possible to simulate reactions involving complex atomic rearrangements and very large energy barriers in runs of a few picoseconds.

A flexible nudged elastic band program for optimization of minimum energy pathways using ab initio electronic structure methods

A driver program for carrying out nudged elastic band optimizations of minimum energy reaction pathways is described. This approach allows for the determination of minimum energy pathways using only energies and gradient information. The driver code has been interfaced with the GAUSSIAN 98 program. Applications to two isomerization reactions and to a cluster model for H(2) desorption from the Si(100)-2 x 1 surface are presented.

Multiple time scale simulations of metal crystal growth reveal the importance of multiatom surface processes

A method for extending atomistic computer simulations of solids beyond the nanosecond time scale was used to simulate metal crystal growth on the time scale of laboratory experiments. Transitions involving concerted motion of multiple atoms on the crystal surface are found to lead to remarkably smooth growth of pure Al(100). Cu(100) is found to grow with a rougher surface, consistent with experiments. Not only is the activation energy of the multiatom Al processes surprisingly low, but vibrational entropy also favors processes where many atoms are displaced.

A first-principles study of methanol decomposition on Pt(111)

A periodic, self-consistent, Density Functional Theory study of methanol decomposition on Pt(111) is presented. The thermochemistry and activation energy barriers for all the elementary steps, starting with O[bond]H scission and proceeding via sequential hydrogen abstraction from the resulting methoxy intermediate, are presented here. The minimum energy path is represented by a one-dimensional potential energy surface connecting methanol with its final decomposition products, CO and hydrogen gas. It is found that the rate-limiting step for this decomposition pathway is the abstraction of hydroxyl hydrogen from methanol. CO is clearly identified as a strong thermodynamic sink in the reaction pathway while the methoxy, formaldehyde, and formyl intermediates are found to have low barriers to decomposition, leading to very short lifetimes for these intermediates. Stable intermediates and transition states are found to obey gas-phase coordination and bond order rules on the Pt(111) surface.

Finding transition paths and rate coefficients through accelerated Langevin dynamics

We present a technique to resolve the rare event problem for a Langevin equation describing a system with thermally activated transitions. A transition event within a given time interval (0,t(f)) can be described by a transition path that has an activation part during (0,t(M)) and a deactivation part during (t(M),t(f))(0<t(M)<t(f)). The activation path is governed by a Langevin equation with negative friction while the deactivation path by the standard Langevin equation with positive friction. Each transition path carries a given statistical weight from which rate constants and related physical quantities can be obtained as averages over all possible paths. We demonstrate how this technique can be used to calculate activation rates of a particle in a two dimensional potential for a wide range of temperatures where standard molecular dynamics techniques are inefficient.

Closing the gap between experiment and theory: crystal growth by temperature accelerated dynamics

We present atomistic simulations of crystal growth where realistic experimental deposition rates are reproduced, without needing any a priori information on the relevant diffusion processes. Using the temperature accelerated dynamics method, we simulate the deposition of 4 monolayers (ML) of Ag/Ag(100) at the rate of 0.075 ML/s, thus obtaining a boost of several orders of magnitude with respect to ordinary molecular dynamics. In the temperature range analyzed (0-70 K), steering and activated mechanisms compete in determining the surface roughness.

Theoretical calculations of dissociative adsorption of CH4 on an Ir(111) surface

The activation energy for chemisorption of CH(4) at an Ir(111) surface is determined using density functional theory combined with an estimate of the long range dispersion interaction. The results are found to be in good agreement with published results of bulb and beam experiments analyzed with a precursor model. A surprisingly large surface relaxation is found where an Ir surface atom is displaced outwards by as much as 0.6 A. A strongly bound molecular state at kinks and adatoms involving eta(2)-H,H bonding was also found.

Theoretical studies of atomic-scale processes relevant to crystal growth

The study of adsorption, diffusion, island formation, and interlayer transport of atoms on a growing surface has been an active field in the past decade, because of both experimental and theoretical advances. Experiments can give detailed images of patterns formed on growing surfaces. An important challenge to the theoretical studies is the identification of dynamical processes controlling the pattern formation and overall surface morphology. This can be achieved by accurate modeling of the atomic interactions, a thorough search for active atomic-scale processes, and simulation of the growth on the experimental timescale to allow for detailed comparison with the experimental measurements. An overview of some of the theoretical methodology used in these studies and results obtained for one of the most extensively studied systems, Pt(111), is given here. A remarkable richness of phenomena has emerged from these studies, where apparently small effects can shift the balance between competing molecular processes and thereby change the morphology of a growing surface. The article concludes with a discussion of possible future directions in this research area.