A method for finding the ridge between saddle points applied to rare event rate estimates

A method is presented for finding the ridge between first order saddle points on a multidimensional surface. For atomic scale systems, such saddle points on the energy surface correspond to atomic rearrangement mechanisms. Information about the ridge can be used to test the validity of the harmonic approximation to transition state theory, in particular to verify that second order saddle points--maxima along the ridge--are high enough compared to the first order saddle points. New minima along the ridge can also be identified during the path optimisation, thereby revealing additional transition mechanisms. The method is based on a string of discretisation points along a path between the first order saddle points and using an iterative optimisation which requires only the force acting on the atoms. At each iteration during the optimisation, the force is inverted along an unstable eigenmode perpendicular to the path. The method is applied to Al adatom diffusion on the Al(100) surface to find the ridge between 2-, 3- and 4-atom concerted displacements and hop mechanisms. A correction to the harmonic approximation of transition state theory was estimated by direct evaluation of the configuration integral along the ridge.     

Comparison of methods for finding saddle points without knowledge of the final states

Within the harmonic approximation to transition state theory, the biggest challenge involved in finding the mechanism or rate of transitions is the location of the relevant saddle points on the multidimensional potential energy surface. The saddle point search is particularly challenging when the final state of the transition is not specified. In this article we report on a comparison of several methods for locating saddle points under these conditions and compare, in particular, the well-established rational function optimization (RFO) methods using either exact or approximate Hessians with the more recently proposed minimum mode following methods where only the minimum eigenvalue mode is found, either by the dimer or the Lanczos method. A test problem involving transitions in a seven-atom Pt island on a Pt(111) surface using a simple Morse pairwise potential function is used and the number of degrees of freedom varied by varying the number of movable atoms. In the full system, 175 atoms can move so 525 degrees of freedom need to be optimized to find the saddle points. For testing purposes, we have also restricted the number of movable atoms to 7 and 1. Our results indicate that if attempting to make a map of all relevant saddle points for a large system (as would be necessary when simulating the long time scale evolution of a thermal system) the minimum mode following methods are preferred. The minimum mode following methods are also more efficient when searching for the lowest saddle points in a large system, and if the force can be obtained cheaply. However, if only the lowest saddle points are sought and the calculation of the force is expensive but a good approximation for the Hessian at the starting position of the search can be obtained at low cost, then the RFO approaches employing an approximate Hessian represent the preferred choice. For small and medium sized systems where the force is expensive to calculate, the RFO approaches employing an approximate Hessian is also the more efficient, but when the force and Hessian can be obtained cheaply and only the lowest saddle points are sought the RFO approach using an exact Hessian is the better choice. These conclusions have been reached based on a comparison of the total computational effort needed to find the saddle points and the number of saddle points found for each of the methods. The RFO methods do not perform very well with respect to the latter aspect, but starting the searches further away from the initial minimum or using the hybrid RFO version presented here improves this behavior considerably in most cases.     

Simplified and improved string method for computing the minimum energy paths in barrier-crossing events

We present a simplified and improved version of the string method, originally proposed by E et al. [Phys. Rev. B 66, 052301 (2002)] for identifying the minimum energy paths in barrier-crossing events. In this new version, the step of projecting the potential force to the direction normal to the string is eliminated and the full potential force is used in the evolution of the string. This not only simplifies the numerical procedure, but also makes the method more stable and accurate. We discuss the algorithmic details of the improved string method, analyze its stability, accuracy and efficiency, and illustrate it via numerical examples. We also show how the string method can be combined with the climbing image technique for the accurate calculation of saddle points and we present another algorithm for the accurate calculation of the unstable directions at the saddle points.     

Force reversed method for locating transition states

To identify the transition state accurately and efficiently on a high-dimensional potential energy surface is one of the most important topics in kinetic studies on chemical reactions. We present here an algorithm to search the transition state by so-called force reversed method, which only requires a rough reaction direction instead of knowing the initial state and final state. Compared to the nudged elastic band method and the dimer method that require multiple images, the present algorithm with only single image required saves significantly the computational cost. The algorithm was implemented in the first-principle periodic total energy calculation package and applied successfully to several prototype surface processes such as the adsorbate diffusion and dissociation on metal surfaces. The results indicate that the force reversed method is efficient, robust to identify the transition state of various surface processes.     

Heating and flooding: a unified approach for rapid generation of free energy surfaces

We propose a general framework for the efficient sampling of conformational equilibria in complex systems and the generation of associated free energy hypersurfaces in terms of a set of collective variables. The method is a strategic synthesis of the adiabatic free energy dynamics approach, previously introduced by us and others, and existing schemes using Gaussian-based adaptive bias potentials to disfavor previously visited regions. In addition, we suggest sampling the thermodynamic force instead of the probability density to reconstruct the free energy hypersurface. All these elements are combined into a robust extended phase-space formalism that can be easily incorporated into existing molecular dynamics packages. The unified scheme is shown to outperform both metadynamics and adiabatic free energy dynamics in generating two-dimensional free energy surfaces for several example cases including the alanine dipeptide in the gas and aqueous phases and the met-enkephalin oligopeptide. In addition, the method can efficiently generate higher dimensional free energy landscapes, which we demonstrate by calculating a four-dimensional surface in the Ramachandran angles of the gas-phase alanine tripeptide.     

Single-sweep methods for free energy calculations

A simple, efficient, and accurate method is proposed to map multidimensional free energy landscapes. The method combines the temperature-accelerated molecular dynamics (TAMD) proposed in [L. Maragliano and E. Vanden-Eijnden, Chem. Phys. Lett. 426, 168 (2006)] with a variational reconstruction method using radial-basis functions for the representation of the free energy. TAMD is used to rapidly sweep through the important regions of the free energy landscape and to compute the gradient of the free energy locally at points in these regions. The variational method is then used to reconstruct the free energy globally from the mean force at these points. The algorithmic aspects of the single-sweep method are explained in detail, and the method is tested on simple examples and used to compute the free energy of the solvated alanine dipeptide in two and four dihedral angles.     

Atomistic simulations of rare events using gentlest ascent dynamics

The dynamics of complex systems often involve thermally activated barrier crossing events that allow these systems to move from one basin of attraction on the high dimensional energy surface to another. Such events are ubiquitous, but challenging to simulate using conventional simulation tools, such as molecular dynamics. Recently, E and Zhou [Nonlinearity 24(6), 1831 (2011)] proposed a set of dynamic equations, the gentlest ascent dynamics (GAD), to describe the escape of a system from a basin of attraction and proved that solutions of GAD converge to index-1 saddle points of the underlying energy. In this paper, we extend GAD to enable finite temperature simulations in which the system hops between different saddle points on the energy surface. An effective strategy to use GAD to sample an ensemble of low barrier saddle points located in the vicinity of a locally stable configuration on the high dimensional energy surface is proposed. The utility of the method is demonstrated by studying the low barrier saddle points associated with point defect activity on a surface. This is done for two representative systems, namely, (a) a surface vacancy and ad-atom pair and (b) a heptamer island on the (111) surface of copper.     

Periodic orbits in biological molecules: phase space structures and selectivity in alanine dipeptide

Small and large molecules may localize their energy in specific bonds or generally in vibrational modes for extended periods of time, an effect which may have dramatic consequences in reaction dynamics. Periodic orbits offer the means to identify phase space regions with localized motions. The author demonstrate that techniques to locate periodic orbits developed for small molecules can be applied to large molecules such as alanine dipeptide. The widely used empirical force fields are employed and principal families of periodic orbits associated with local-type motions and emanated from the lowest energy minima and saddle points are investigated. Continuation of these families at high energies unravels the stable and unstable regions of phase space as well as elementary bifurcations such as saddle nodes.     

分子动力学模拟多巴在自组装膜上的黏附性

为了更好地理解贻贝在表面的黏附机理,实现水下胶黏,采用分子动力学方法研究了多巴在自组装膜上的黏附性:采用伞形取样和加权柱状图分析方法计算了多巴在不同自组装膜表面的黏附自由能,使用拉伸分子动力学模拟研究了多巴在不同自组装膜表面上黏附后的脱附力.结果表明,多巴在带负电的羧基自组装膜上的黏附能比在带正电的氨基自组装膜上的大,多巴更容易黏附到带负电表面;多巴在带电表面的黏附能比未带电表面的黏附能更强,表明在带电表面黏附更稳定.进一步分析了多巴在不同表面的取向分布,发现多巴与不同表面相互作用的方式不同:与疏水表面主要通过苯环相互作用;与亲水表面主要通过羟基相互作用;与负电表面主要通过氨基相互作用;与正电表面主要通过羧基相互作用.通过模拟比较了多巴在不同自组装膜上的脱附力,发现多巴在带电表面的脱附力比在未带电表面的大,与黏附能的趋势一致.对比4种非带电表面的脱附力,发现多巴在疏水性甲基自组装膜表面的脱附力最大,黏附更稳定,随着表面疏水性的增加,脱附力增大,黏附稳定性增强.本工作可为研发新型水下胶黏剂提供理论指导.     

Prediction of the orientations of adsorbed protein using an empirical energy function with implicit solvation

When simulating protein adsorption behavior, decisions must first be made regarding how the protein should be oriented on the surface. To address this problem, we have developed a molecular simulation program that combines an empirical adsorption free energy function with an efficient configurational search method to calculate orientation-dependent adsorption free energies between proteins and functionalized surfaces. The configuration space is searched systematically using a quaternion rotation technique, and the adsorption free energy is evaluated using an empirical energy function with an efficient grid-based calculational method. In this paper, the developed method is applied to analyze the preferred orientations of a model protein, lysozyme, on various functionalized alkanethiol self-assembled monolayer (SAM) surfaces by the generation of contour graphs that relate adsorption free energy to adsorbed orientation, and the results are compared with experimental observations. As anticipated, the adsorbed orientation of lysozyme is predicted to be dependent on the discrete organization of the functional groups presented by the surface. Lysozyme, which is a positively charged protein, is predicted to adsorb on its 'side' on both hydrophobic and negatively charged surfaces. On surfaces with discrete positively charged sites, attractive interaction energies can also be obtained due to the presence of discrete local negative charges present on the lysozyme surface. In this case, 'end-on' orientations are preferred. Additionally, SAM surface models with mixed functionality suggest that the interactions between lysozyme and surfaces could be greatly enhanced if individual surface functional groups are able to access the catalytic cleft region of lysozyme, similar to ligand-receptor interactions. The contour graphs generated by this method can be used to identify low-energy orientations that can then be used as starting points for further simulations to investigate conformational changes induced in protein structure following initial adsorption.     

Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules

Many interesting dynamic properties of biological molecules cannot be simulated directly using molecular dynamics because of nanosecond time scale limitations. These systems are trapped in potential energy minima with high free energy barriers for large numbers of computational steps. The dynamic evolution of many molecular systems occurs through a series of rare events as the system moves from one potential energy basin to another. Therefore, we have proposed a robust bias potential function that can be used in an efficient accelerated molecular dynamics approach to simulate the transition of high energy barriers without any advance knowledge of the location of either the potential energy wells or saddle points. In this method, the potential energy landscape is altered by adding a bias potential to the true potential such that the escape rates from potential wells are enhanced, which accelerates and extends the time scale in molecular dynamics simulations. Our definition of the bias potential echoes the underlying shape of the potential energy landscape on the modified surface, thus allowing for the potential energy minima to be well defined, and hence properly sampled during the simulation. We have shown that our approach, which can be extended to biomolecules, samples the conformational space more efficiently than normal molecular dynamics simulations, and converges to the correct canonical distribution.     

Free energy perturbation study of water dimer dissociation kinetics

An efficient approach is described for using accurate ab initio calculations to determine the rates of elementary condensation and evaporation processes that lead to nucleation of aqueous aerosols. The feasibility of the method is demonstrated in an application to evaporation rates of water dimer at 230 K. The method, known as ABC-FEP (ab initio/classical free energy perturbation), begins with a calculation of the potential of mean force for the dissociation (evaporation) of small water clusters using a molecular dynamics (MD) simulation with a model potential. The free energy perturbation is used to calculate how changing from the model potential to a potential calculated from ab initio methods would alter the potential of mean force. The difference in free energy is the Boltzmann-weighted average of the difference between the ab initio and classical potential energies, with the average taken over a sample of configurations from the MD simulation. In principle, the method does not require a highly accurate model potential, though more accurate potentials require fewer configurations to achieve a small sampling error in the free energy perturbation step. To test the feasibility of obtaining accurate potentials of mean force from ab initio calculations at a modest number of configurations, the free energy perturbation method has been used to correct the errors when some standard models for bulk water (SPC, TIP4P, and TIP4PFQ) are applied to water dimer. To allow a thorough exploration of sampling issues, a highly accurate fit to results of accurate ab initio calculations, known as SAPT-5s, as been used a proxy for the ab initio calculations. It is shown that accurate values for a point on the potential of mean force can be obtained from any of the water models using ab initio calculations at only 50 configurations. Thus, this method allows accurate simulations of small clusters without the need to develop water models specifically for clusters.     

Multiple pathways in conformational transitions of the alanine dipeptide: an application of dynamic importance sampling

We present multiple dynamic transition pathways on the two-dimensional dihedral plane between conformational states of the alanine dipeptide. The method used in this study is dynamic importance sampling (DIMS). To perform DIMS, unbiased molecular dynamic simulations are used to generate equilibrium ensembles for the alanine dipeptide within different states. Free energy surfaces on the dihedral plane are calculated from the equilibrium simulations, and four energy minima defined from the surface are used as the starting and ending points for DIMS dynamics. The DIMS method represents an important step towards finding multiple transition pathways within complex biomolecular systems.     

Simple physics-based analytical formulas for the potentials of mean force for the interaction of amino acid side chains in water. 1. Approximate expression for the free energy of hydrophobic association based on a Gaussian-overlap model

A physics-based model is proposed to derive approximate analytical expressions for the cavity component of the free energy of hydrophobic association of spherical and spheroidal solutes in water. The model is based on the difference between the number and context of the water molecules in the hydration sphere of a hydrophobic dimer and of two isolated hydrophobic solutes. It is assumed that the water molecules touching the convex part of the molecular surface of the dimer and those in the hydration spheres of the monomers contribute equally to the free energy of solvation, and those touching the saddle part of the molecular surface of the dimer result in a more pronounced increase in free energy because of their more restricted mobility (entropy loss) and fewer favorable electrostatic interactions with other water molecules. The density of water in the hydration sphere around a single solute particle is approximated by the derivative of a Gaussian centered on the solute molecule with respect to its standard deviation. On the basis of this approximation, the number of water molecules in different parts of the hydration sphere of the dimer is expressed in terms of the first and the second mixed derivatives of the two Gaussians centered on the first and second solute molecules, respectively, with respect to the standard deviations of these Gaussians, and plausible analytical expressions for the cavity component of the hydrophobic-association energy of spherical and spheroidal solutes are introduced. As opposed to earlier hydration-shell models, our expressions reproduce the desolvation maxima in the potentials of mean force of pairs of nonpolar solutes in water, and their advantage over the models based on molecular-surface area is that they have continuous gradients in the coordinates of solute centers.     

Lowest singlet and triplet potential energy surfaces of S2N2

Forty four stationary points have been located on the lowest singlet and triplet potential energy surfaces of S(2)N(2). Ten minima and ten saddle points on the lowest singlet surface and eleven minima and thirteen saddle points on the lowest triplet surface were found. All saddle points were connected to minima or lower-order saddle points by following the intrinsic reaction coordinate. Renner-Teller effects in the linear isomers were studied by examining their bending curves. The S(2)N(2) polymerization mechanism was investigated by first locating the transition state corresponding to ring opening and then considering all species connected to it that are close in energy. The commonly accepted mechanism is problematic due to the number of species that would lead to dissociation to SN + SN. Other possible isomers that are consistent with the experimental evidence but do not connect to SN radicals in the dissociation limit were examined. A mechanism of polymerization to (SN)(x)() is proposed that involves excitation of the square planar singlet molecule to the triplet surface. The triplet species then undergoes a puckering, and polymerization occurs in a direction approximately perpendicular to the S(2)N(2) plane. Consideration of the predicted vibrational frequencies suggests the structure of the second isomer of S(2)N(2). This isomer has a trans-NSSN structure with a long SS bond. The energetics of trans-NSSN are consistent with the observed temperature effects in the dimerization of SN. Analysis of the bending curves of linear NSSN and NSNS indicates that trans-NSSN is the only isomer which has a small yet significant barrier to that dimerization.     

A new form of kinetic isotope effect. Dynamic effects on isotopic selectivity and regioselectivity

The intramolecular H/D kinetic isotope effect in the ene reaction of singlet oxygen with tetramethylethylene is studied using quasiclassical direct dynamics calculations on a B3LYP/6-31G* potential energy surface. Starting from the area of the energy surface around a valley-ridge inflection point, random trajectories lead to predominantly H abstraction over D abstraction, despite the symmetry of the surface and the absence of a barrier to either reaction. This demonstrates a new form of kinetic isotope effect, unrelated to the usual effect of zero-point energies on barriers. Dynamics calculations on the reaction of cis-2-pentene predict the experimentally observed mixture of regioisomeric products, while the minimum-energy path leads to only one product. For energy surfaces containing two adjacent saddle points, dynamics effects are important for understanding both product and isotopic selectivity, and this should be considered in the interpretation of experimental results.     

Calibration of force-field dependency in free energy landscapes of peptide conformations by quantum chemical calculations

The free energy landscapes of peptide conformations were calibrated by ab initio quantum chemical calculations, after the enhanced conformational diversity search using the multicanonical molecular dynamics simulations. Three different potentials of mean force for an isolated dipeptide were individually obtained by the multicanonical molecular dynamics simulations using the conventional force fields, AMBER parm94, AMBER parm96, and CHARMm22. Each potential of mean force was then calibrated based upon the umbrella sampling algorithm from the adiabatic energy map that was calculated separately by the ab initio molecular orbital method, and all of the calibrated potentials of mean force coincided well. The calibration method was also applied to the simulations of a peptide dimer in explicit water models, and it was shown that the calibrated free energy landscapes did not depend on the force field used in the classical simulations, as far as the conformational space was sampled well. The current calibration method fuses the classical free energy calculation with the quantum chemical calculation, and it should generally make simulations for biomolecular systems much more reliable when combining with enhanced conformational sampling.     

Optimizing transition states via kernel-based machine learning

We present a method for optimizing transition state theory dividing surfaces with support vector machines. The resulting dividing surfaces require no a priori information or intuition about reaction mechanisms. To generate optimal dividing surfaces, we apply a cycle of machine-learning and refinement of the surface by molecular dynamics sampling. We demonstrate that the machine-learned surfaces contain the relevant low-energy saddle points. The mechanisms of reactions may be extracted from the machine-learned surfaces in order to identify unexpected chemically relevant processes. Furthermore, we show that the machine-learned surfaces significantly increase the transmission coefficient for an adatom exchange involving many coupled degrees of freedom on a (100) surface when compared to a distance-based dividing surface.     

Direct dynamics simulations of O(3P) + HCl at hyperthermal collision energies

The dynamics of the O(3P) + HCl reaction at hyperthermal collision energies were investigated using the quasiclassical trajectory method. Stationary points on the OClH 3A" and 3A' potential energy surfaces (PESs) were also examined. The lowest transition state leading to OCl + H on the 3A" surface is 2.26 eV above the reagents at the CCSD(T)/cc-pVTZ level of theory. This saddle point is bent and product-like. Direct dynamics calculations at the MP2/cc-pVTZ level of theory were used to investigate the excitation functions for OH + Cl, OCl + H, and O + H + Cl formation. OCl is formed mainly from small-impact-parameter collisions, and the OCl + H excitation function peaks around 5 eV, where it is similar in magnitude to the OH + Cl excitation function. The shape of the OCl + H excitation function is discussed, and features are identified that should be general to hyperthermal collision dynamics.