Computer optimization of small atomic clusters

We have successfully identified stable configurations of both rare-gas and NaCl clusters with a new optimization procedure. An initial cluster configuration is prepared in a so-called shoot-and-stay process. Its total energy is then minimized with respect to the atomic coordinates. To prevent the system from being locked in local minima, the step size of each move is chosen as the width of the energy well at a higher level. As the system evolves, the global minimum is contained in the volume bounded by the decreasing value of step sizes. We have also carried out the optimization of NaCl clusters by the simulated annealing technique, for comparison. The results show that for such heterogeneous systems, the latter method cannot always find the global minimum, because of large energy gaps between different catchment regions in phase space.     

Equilibrium thermodynamics from basin-sampling

We present a "basin-sampling" approach for calculation of the potential energy density of states for classical statistical models. It combines a Wang-Landau-type uniform sampling of local minima and a novel approach for approximating the relative contributions from local minima in terms of the volumes of basins of attraction. We have employed basin-sampling to study phase changes in atomic clusters modeled by the Lennard-Jones potential and for ionic clusters. The approach proves to be efficient for systems involving broken ergodicity and has allowed us to calculate converged heat capacity curves for systems that could previously only be treated using the harmonic superposition approximation. Benchmarks are also provided by comparison with parallel tempering and Wang-Landau simulations, where these proved feasible.     

Rotational symmetry of the molecular potential energy in the Cartesian coordinates

We consider the molecular Born-Oppenheimer potential energy as a function of atomic Cartesian coordinates and discuss the non-stationary Hessian properties arising due to rotational symmetry. A connection with the extended Hessian theory is included. New applications of Cartesian representation for examining and correcting raw numerical Hessian data and a simple formulation of harmonic vibrational analysis of partially optimized systems are proposed. Exemplary calculations for the porphyrin molecule with an internal proton transfer are reported. We also develop the normal transformation method to incorporate the rotational symmetry into the approximate analytical potentials, which are parametrized in the Cartesian coordinates. The transformation converts the coordinates from the space fixed frame to the frame which translates and rotates with the molecule and is determined by the Eckart conditions. New simple analytical formulas for the first and second derivatives of the transformed potential are derived. This fast method can be used to calculate the potential and its derivatives in the simulations of chemical reaction dynamics in the space fixed Cartesian frame without the need to constrain the molecular rotation or to define the local non-redundant internal coordinates.     

Atomic mobility in strained glassy polymers: The role of fold catastrophes on the potential energy surface

Deformation is known to enhance the atomic mobility in disordered systems such as polymer materials. To elucidate the origin of this phenomenon, we carry out two types of simulations: molecular dynamics (MD) simulations, which determine the atomic trajectories at finite temperature, and quasi-static simulations, which determine the atomic trajectories in the limit of zero temperature (and in the limit of zero shear rate). The quasi-static simulations show discontinuous changes in properties, such as system energy and atomic mobility. We use a normal mode analysis to show that these discontinuous changes arise from fold catastrophes of the potential energy landscape, in which energy minima flatten out and the heights of energy barriers decrease to zero; this was demonstrated by normal mode frequencies following a power law with an exponent of 0.5 as the discontinuous change is approached. After the fold catastrophe, the system relaxes to a different energy minimum, giving rise to atomic displacements. These fold catastrophes are the only mechanism for diffusive atomic displacements in the quasi-static simulations, where there is no thermal energy. We compared the mean-squared displacements as a function of strain from the quasi-static simulations to those from MD simulations (which do include thermal effects)—the similarity of the values of the mean-squared displacements in these two types of simulations demonstrates that the fold catastrophes underlie the enhanced dynamics in strained polymer systems even at finite temperature. © 2012 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys, 2012     

Toward absolute density of states calculations for proteins

The density of states (DOS), which gives the number of conformations with a particular energy E, is a prerequisite in computing most thermodynamic quantities and in elucidating important biological processes such as the mechanism of protein folding. However, current methods for computing DOS of large systems such as proteins generally yield only the ratios of microstate counts for different energies, which could yield absolute conformation counts if the total number of conformations in phase space is known, thus motivating this work. Here, the total number of energy minima of 50-mer polyalanine, whose size corresponds to naturally occurring small proteins, was estimated under an all-atom potential energy function based on the cumulative distribution function (CDF) of conformational differences to be approximately 10(38). This estimate can place any DOS function, such as the Gaussian DOS distribution in the random energy model, on an absolute scale. Comparing the absolute conformational counts from a Gaussian DOS function with those from the CDF derived from quenched molecular dynamics ensembles shows that the former are far greater than the latter, indicating far fewer low-energy minima actually exist. In addition to showing how CDF and relative DOS calculations can yield absolute DOS for a discrete system, we also show how they can yield absolute DOS for continuous variable systems to a specified atomic variance. In the context of protein folding, knowing this phase-space "volume" of conformations in a DOS function, as well as characteristic transition times, constrains the set of possible folding mechanisms.     

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.     

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.     

Accelerating molecular simulations by reversible mapping between local minima

A new framework is presented for performing Monte Carlo simulations of condensed matter based on a recently developed bijective mapping between local energy minima. The framework is used to implement a range of new multiparticle Monte Carlo moves, which are investigated by simulating atomic Lennard-Jones fluids in the canonical and grand canonical ensembles. Important aspects of the simulation protocol and their effect on performance are analyzed in detail. Using the mapping accelerates the simulations by many orders of magnitude when compared to the equivalent moves without the mapping, and leads to particularly efficient configurational sampling at low temperatures and high densities. The method appears to be suitable for adapting to quantitative simulations of more complex molecular systems over long effective time scales.     

Potential energy and free energy landscapes

Familiar concepts for small molecules may require reinterpretation for larger systems. For example, rearrangements between geometrical isomers are usually considered in terms of transitions between the corresponding local minima on the underlying potential energy surface, V. However, transitions between bulk phases such as solid and liquid, or between the denatured and native states of a protein, are normally addressed in terms of free energy minima. To reestablish a connection with the potential energy surface we must think in terms of representative samples of local minima of V, from which a free energy surface is projected by averaging over most of the coordinates. The present contribution outlines how this connection can be developed into a tool for quantitative calculations. In particular, stepping between the local minima of V provides powerful methods for locating the global potential energy minimum, and for calculating global thermodynamic properties. When the transition states that link local minima are also sampled we can exploit statistical rate theory to obtain insight into global dynamics and rare events. Visualizing the potential energy landscape helps to explain how the network of local minima and transition states determines properties such as heat capacity features, which signify transitions between free energy minima. The organization of the landscape also reveals how certain systems can reliably locate particular structures on the experimental time scale from among an exponentially large number of local minima. Such directed searches not only enable proteins to overcome Levinthal's paradox but may also underlie the formation of "magic numbers" in molecular beams, the self-assembly of macromolecular structures, and crystallization.     

Conformational analysis of a stereochemically complete set of cis-enediol peptide analogues

A conformational analysis of a stereochemically complete set of peptide analogues based on a cis-enediol unit is presented. The cis-enediol unit, which can replace a two or a three amino acid segment of a peptide, contains two "side chains", four asymmetrical carbon atoms, and six free dihedral angles. To determine the accessible conformational space, the molecules are divided into three fragments, each containing two free dihedral angles. The energy surfaces are computed for all dihedral angle values, and the possible conformations of the cis-enediol unit analogues are built using all combinations of the surface minima. Such a "build-up" procedure, which is very fast, is able to reproduce 75% of the minima obtained from a full dihedral angle exploration of the conformational space. The cis-enediol unit minima are compared with the corresponding di- and tripeptide minima; all peptide minima can be closely matched by a cis-enediol unit minimum of low energy (less than 2.2 kcal/mol above the lowest energy conformer). However, there are low energy minima of the cis-enediol unit that have no corresponding minima in peptides. The results are shown to depend strongly on the chirality of the analogues. The ability of each of the stereoisomers to mimic natural peptides, evaluated by the present approach, is correlated with its experimental activity in a renin inhibition assay.     

Generating conformer ensembles using a multiobjective genetic algorithm

The task of generating a nonredundant set of low-energy conformations for small molecules is of fundamental importance for many molecular modeling and drug-design methodologies. Several approaches to conformer generation have been published. Exhaustive searches suffer from the exponential growth of the search space with increasing degrees of conformational freedom (number of rotatable bonds). Stochastic algorithms do not suffer as much from the exponential increase of search space and provide a good coverage of the energy minima. Here, the use of a multiobjective genetic algorithm in the generation of conformer ensembles is investigated. Distance geometry is used to generate an initial conformer, which is then subject to geometric modifications encoded by the individuals of the genetic algorithm. The geometric modifications apply to torsion angles about rotatable bonds, stereochemistry of double bonds and tetrahedral chiral centers, and ring conformations. The geometric diversity of the evolving conformer ensemble is preserved by a fitness-sharing mechanism based on the root-mean-square distance of the atomic coordinates. Molecular symmetry is taken into account in the distance calculation. The geometric modifications introduce strain into the structures. The strain is relaxed using an MMFF94-like force field in a postprocessing step that also removes conformational duplicates and structures whose strain energy remains above a predefined window from the minimum energy value found in the set. The implementation, called Balloon, is available free of charge on the Internet ( http://www.abo.fi/~mivainio/balloon/).     

Conformational distribution of heptaalanine: Analysis using a new Monte Carlo chain growth method

A distribution of conformations of heptaalanine is obtained using a new Monte Carlo (MC) method that grows the chain atom by atom. Resulting configurations are energy minimized and a detailed analysis is performed of the minimum-energy conformers using a method of classification that partitions ?ψ space. The MC-generated configurations are compared with those generated from high-temperature molecular dynamics (MD) simulations. It is found that the new method generates a wide distribution of low-energy conformers at least 10 times more quickly than the MD. An analysis of the generated energy minima demonstrates that they can be divided into clusters in the space defined by the five pairs of ?—ψ angles of the inner residues. The space occupied by the minima populating each cluster is restricted. © 1992 by John Wiley & Sons, Inc.     

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.     

Global optimization of ionic Mg(n)F(2n) (n=1-30) clusters

The global optimization basin-hopping (BH) method has been used to locate the global minima (GM) of Mg(n)F(2n) (n=1-30) clusters using a Born-Mayer-type potential. Some of the GM were particularly difficult to find, requiring more than 1.5 x 10(4) BH steps. We have found that both the binding energy per MgF2 unit and the effective volume of the GM isomers increase almost linearly with n, and that cluster symmetry decreases with cluster size. The data derived from the BH runs reveal a growing density of local minima just above the GM as n increases. Despite this, the attraction basin around each GM is relatively large, since after all their atomic coordinates are randomly displaced by values as high as 2.0 bohrs, the perturbed structures, upon reoptimization, relax back to the GM in more than 50% of the cases (except for n=10 and 11). The relative stabilities derived from energy second differences suggest that n=8,10,13,15, and 20 are probably the magic numbers for these systems. Mass spectrum experiments would be very useful to clarify this issue.     

Overcoming entropic barrier with coupled sampling at dual resolutions

An enhanced sampling method is proposed for ab initio protein folding simulations. The new method couples a high-resolution model for accuracy and a low-resolution model for efficiency. It aims to overcome the entropic barrier found in the exponentially large protein conformational space when a high-resolution model, such as an all-atom molecular mechanics force field, is used. The proposed method is designed to satisfy the detailed balance condition so that the Boltzmann distribution can be generated in all sampling trajectories in both high and low resolutions. The method was tested on model analytical energy functions and ab initio folding simulations of a beta-hairpin peptide. It was found to be more efficient than replica-exchange method that is used as its building block. Analysis with the analytical energy functions shows that the number of energy calculations required to find global minima and to converge mean potential energies is much fewer with the new method. Ergodic measure shows that the new method explores the conformational space more rapidly. We also studied imperfect low-resolution energy models and found that the introduction of errors in low-resolution models does decrease its sampling efficiency. However, a reasonable increase in efficiency is still observed when the global minima of the low-resolution models are in the vicinity of the global minimum basin of the high-resolution model. Finally, our ab initio folding simulation of the tested peptide shows that the new method is able to fold the peptide in a very short simulation time. The structural distribution generated by the new method at the equilibrium portion of the trajectory resembles that in the equilibrium simulation starting from the crystal structure.     

Minima hopping: an efficient search method for the global minimum of the potential energy surface of complex molecular systems

A method is presented that can find the global minimum of very complex condensed matter systems. It is based on the simple principle of exploring the configurational space as fast as possible and of avoiding revisiting known parts of this space. Even though it is not a genetic algorithm, it is not based on thermodynamics. The efficiency of the method depends strongly on the type of moves that are used to hop into new local minima. Moves that find low-barrier escape-paths out of the present minimum generally lead into low energy minima.     

On the Use of Accelerated Molecular Dynamics to Enhance Configurational Sampling in Ab Initio Simulations

We have implemented the accelerated molecular dynamics approach (Hamelberg, D.; Mongan, J.; McCammon, J. A. J. Chem. Phys. 2004, 120 (24), 11919) in the framework of ab initio MD (AIMD). Using three simple examples, we demonstrate that accelerated AIMD (A-AIMD) can be used to accelerate solvent relaxation in AIMD simulations and facilitate the detection of reaction coordinates: (i) We show, for one cyclohexane molecule in the gas phase, that the method can be used to accelerate the rate of the chair-to-chair interconversion by a factor of ～1 × 10(5), while allowing for the reconstruction of the correct canonical distribution of low-energy states; (ii) We then show, for a water box of 64 H(2)O molecules, that A-AIMD can also be used in the condensed phase to accelerate the sampling of water conformations, without affecting the structural properties of the solvent; and (iii) The method is then used to compute the potential of mean force (PMF) for the dissociation of Na-Cl in water, accelerating the convergence by a factor of ～3-4 compared to conventional AIMD simulations.(2) These results suggest that A-AIMD is a useful addition to existing methods for enhanced conformational and phase-space sampling in solution. While the method does not make the use of collective variables superfluous, it also does not require the user to define a set of collective variables that can capture all the low-energy minima on the potential energy surface. This property may prove very useful when dealing with highly complex multidimensional systems that require a quantum mechanical treatment.     

Estimating kinetic rates from accelerated molecular dynamics simulations: alanine dipeptide in explicit solvent as a case study

Molecular dynamics (MD) simulation is the standard computational technique used to obtain information on the time evolution of the conformations of proteins and many other molecular systems. However, for most biological systems of interest, the time scale for slow conformational transitions is still inaccessible to standard MD simulations. Several sampling methods have been proposed to address this issue, including the accelerated molecular dynamics method. In this work, we study the extent of sampling of the phi/psi space of alanine dipeptide in explicit water using accelerated molecular dynamics and present a framework to recover the correct kinetic rate constant for the helix to beta-strand transition. We show that the accelerated MD can drastically enhance the sampling of the phi/psi conformational phase space when compared to normal MD. In addition, the free energy density plots of the phi/psi space show that all minima regions are accurately sampled and the canonical distribution is recovered. Moreover, the kinetic rate constant for the helix to beta-strand transition is accurately estimated from these simulations by relating the diffusion coefficient to the local energetic roughness of the energy landscape. Surprisingly, even for such a low barrier transition, it is difficult to obtain enough transitions to accurately estimate the rate constant when one uses normal MD.     

Equilibrium density of states and thermodynamic properties of a model glass former

This paper presents an analysis of the thermodynamics of a model glass former. We have performed equilibrium sampling of a popular binary Lennard-Jones model, employing parallel tempering Monte Carlo to cover the crystalline, amorphous, and liquid regions of configuration space. Disconnectivity graphs are used to visualize the potential energy landscape in the vicinity of a crystalline geometry and in an amorphous region of configuration space. The crystalline global minimum is separated from the bulk of the minima by a large potential energy gap, leading to broken ergodicity in conventional simulations. Our sampling reveals crystalline global minima that are lower in potential energy than some of the previous candidates. We present equilibrium thermodynamic properties based on parallel tempering simulations, including heat capacities and free energy profiles, which depend explicitly on the crystal structure. We also report equilibrium melting temperatures.     

Metadynamics in essential coordinates: free energy simulation of conformational changes

We propose an approach that combines an extraction of collective motions of a molecular system with a sampling of its free energy surface. A recently introduced method of metadynamics allows exploration of the free energy surface of a molecular system by means of coarse-grained dynamics with flooding of free energy minima. This free energy surface is defined as a function of a set of collective variables (e.g., interatomic distances, angles, torsions, and others). In this study, essential coordinates determined by essential dynamics (principle component analysis) were used as collective variables in metadynamics. First, dynamics of the model system (explicitly solvated alanine dipeptide, Ace-Ala-Nme) was simulated by a classical molecular dynamics simulation. The trajectory (1 ns) was then analyzed by essential dynamics to obtain essential coordinates. The free energy surface as a function of the first and second essential coordinates was then explored by metadynamics. The resulting free energy surface is in agreement with other studies of this system. We propose that a combination of these two methods (metadynamics and essential dynamics) has great potential in studies of conformational changes in peptides and proteins.