A consensus line search algorithm for molecular potential energy functions

Force field based energy minimization of molecular structures is a central task in computational chemistry and biology. Solving this problem usually requires efficient local minimization techniques, i.e., iterative two‐step methods that search first for a descent direction and then try to estimate the step width. The second step, the so called line search, typically uses polynomial interpolation schemes to estimate the next trial step. However, dependent on local properties of the objective function alternative schemes may be more appropriate especially if the objective function shows singularities or exponential behavior. As the choice of the best interpolation scheme cannot be made a priori, we propose a new consensus line search approach that performs several different interpolation schemes at each step and then decides which one is the most reliable at the current position. Although a naive consensus approach would lead to severe performance impacts, our method does not require additional evaluations of the energy function, imposing only negligible computational overhead. Additionally, our method can be easily adapted to the local behavior of other objective functions by incorporating suitable interpolation schemes or omitting non‐fitting schemes. The performance of our consensus line search approach has been evaluated and compared to established standard line search algorithms by minimizing the structures of a large set of molecules using different force fields. The proposed algorithm shows better performance in almost all test cases, i.e., it reduces the number of iterations and function and gradient evaluations, leading to significantly reduced run times. © 2008 Wiley Periodicals, Inc. J Comput Chem, 2009     

Repeated‐annealing sampling combined with multicanonical algorithm for conformational sampling of bio‐molecules

A novel conformational sampling method (repeated‐annealing sampling method) is proposed to execute an efficient conformational sampling at a reasonable computational cost. In the method, a molecular dynamics simulation is done with repeating an elemental process. An elemental process consists of four subprocesses: high‐temperature run, annealing, room‐temperature run, and fast heating. The sampling is done automatically according to a temperature‐control schedule. The room‐temperature run is treated with the multicanonical algorithm, and the other subprocesses are done with the conventional molecular dynamics algorithm. The method, differing from the generalized ensemble methods recently developed, is not warrantable to give the canonical ensemble because of the nonphysical process in the annealing. However, we observed that the slower the annealing and the longer the high‐temperature run, the closer the sampled conformations to those of the canonical ensemble. A test was performed with tri‐N‐acetyl‐D ‐glucosamine in vacuo, and the results were compared with those from the conventional multicanonical simulation. Not only the reweighted canonical distribution function but also the energy landscape were in good agreement with those from the conventional multicanonical simulation. The potential of mean force also showed a fairly good agreement with that from the conventional multicanonical simulation in the room‐temperature region. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1098–1106, 2001     

Molecular structure matching by simulated annealing. III. The incorporation of null correspondences into the matching problem

Summary This paper extends an application of the method of simulated annealing for molecular matching so that the best common subsets of atom positions can be identified. Null correspondences are introduced into the difference distance matrix to enable poorly matched positions to be ignored in minimizing the objective function. The efficiency of the algorithm in finding correct subsets is rigorously tested.     

Time-efficient flexible superposition of medium-sized molecules

We present an efficient algorithm for the structural alignment of medium-sized organic molecules. The algorithm has been developed for applications in 3D QSAR and in receptor modeling. The method assumes one of the molecules, the reference ligand, to be presented in the conformation that it adopts inside the receptor pocket. The second molecule, the test ligand, is considered to be flexible, and is assumed to be given in an arbitrary low-energy conformation. Ligand flexibility is modeled by decomposing the test ligand into molecular fragments, such that ring systems are completely contained in a single fragment. Conformations of fragments and torsional angles of single bonds are taken from a small finite set, which depends on the fragment and bond, respectively. The algorithm superimposes a distinguished base fragment of the test ligand onto a suitable region of the reference ligand and then attaches the remaining fragments of the test ligand in a step-by-step fashion. During this process, a scoring function is optimized that encompasses bonding terms and terms accounting for steric overlap as well as for similarity of chemical properties of both ligands. The algorithm has been implemented in the FLEXS system. To validate the quality of the produced results, we have selected a number of examples for which the mutual superposition of two ligands is experimentally given by the comparison of the binding geometries known from the crystal structures of their corresponding protein–ligand complexes. On more than two-thirds of the test examples the algorithm produces rms deviations of the predicted versus the observed conformation of the test ligand below 1.5 Å. The run time of the algorithm on a single problem instance is a few minutes on a common-day workstation. The overall goal of this research is to drastically reduce run times, while limiting the inaccuracies of the model and the computation to a tolerable level.     

Massive docking of flexible ligands using environmental niches in parallelized genetic algorithms

Virtual screening of large libraries of small compounds requires fast and reliable automatic docking methods. In this article we present a parallel implementation of a genetic algorithm (GA) and the implementation of an enhanced genetic algorithm (EGA) with niching that lead to remarkable speedups compared to the original version AutoDock 3.0. The niching concept is introduced naturally by sharing genetic information between evolutions of subpopulations that run independently, each on one CPU. A unique set of additionally introduced search parameters that control this information flow has been obtained for drug‐like molecules based on the detailed study of three test cases of different complexity. The average docking time for one compound is of 8.6 s using eight R10,000 processors running at 200 MHz in an Origin 2000 computer. Different genetic algorithms with and without local search (LS) have been compared on an equal workload basis showing EGA/LS to be superior over all alternatives because it finds lower energy solutions faster and more often, particularly for high dimensionality problems. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1971–1982, 2001     

An exploration of a novel strategy for superposing several flexible molecules

Summary This paper describes a computational strategy for the superposition of a set of flexible molecules. The combinatorial problems of searching conformational space and molecular matching are reduced drastically by the combined use of simulated annealing methods and cluster analysis. For each molecule, the global minimum of the conformational energy is determined by annealing and the search trajectory is retained in a history file. All the significantly different low-energy conformations are extracted by cluster analysis of data in the history file. Each pair of molecules, in each of their significantly different conformations, is then matched by simulated annealing, using the difference-distance matrix as the objective function. A set of match statistics is then obtained, from which the best match taken from all different conformations can be found. The molecules are then superposed either by reference to a base molecule or by a consensus method. This strategy ensures that as wide a range of conformations as possible is considered, but at the same time that the smallest number of significantly different conformations is used. The method has been tested on a set of six angiotensin II antagonists with between 7–11 rotatable bonds.     

MultiMCS: a fast algorithm for the maximum common substructure problem on multiple molecules

Several efficient correspondence graph-based algorithms for determining the maximum common substructure (MCS) of a pair of molecules have been published in the literature. The extension of the problem to three or more molecules is however nontrivial; heuristics used to increase the efficiency in the two-molecule case are either inapplicable to the many-molecule case or do not provide significant speedups. Our specific algorithmic contribution is two-fold. First, we show how the correspondence graph approach for the two-molecule case can be generalized to obtain an algorithm that is guaranteed to find the optimum connected MCS of multiple molecules, and that runs fast on most families of molecules using a new divide-and-conquer strategy that has hitherto not been reported in this context. Second, we provide a characterization of those compound families for which the algorithm might run slowly, along with a heuristic for speeding up computations on these families. We also extend the above algorithm to a heuristic algorithm to find the disconnected MCS of multiple molecules and to an algorithm for clustering molecules into groups, with each group sharing a substantial MCS. Our methods are flexible in that they provide exquisite control on various matching criteria used to define a common substructure.     

Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function

A novel and robust automated docking method that predicts the bound conformations of flexible ligands to macromolecular targets has been developed and tested, in combination with a new scoring function that estimates the free energy change upon binding. Interestingly, this method applies a Lamarckian model of genetics, in which environmental adaptations of an individual's phenotype are reverse transcribed into its genotype and become heritable traits (sic). We consider three search methods, Monte Carlo simulated annealing, a traditional genetic algorithm, and the Lamarckian genetic algorithm, and compare their performance in dockings of seven protein–ligand test systems having known three-dimensional structure. We show that both the traditional and Lamarckian genetic algorithms can handle ligands with more degrees of freedom than the simulated annealing method used in earlier versions of AUTO DOCK , and that the Lamarckian genetic algorithm is the most efficient, reliable, and successful of the three. The empirical free energy function was calibrated using a set of 30 structurally known protein–ligand complexes with experimentally determined binding constants. Linear regression analysis of the observed binding constants in terms of a wide variety of structure-derived molecular properties was performed. The final model had a residual standard error of 9.11 kJ mol⁻¹ (2.177 kcal mol⁻¹) and was chosen as the new energy function. The new search methods and empirical free energy function are available in AUTO DOCK , version 3.0. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1639–1662, 1998     

Molecular structure matching by simulated annealing. I. A comparison between different cooling schedules

Summary This paper outlines an application of the theory of simulated annealing to molecular matching problems. Three cooling schedules are examined: linear, exponential and dynamic cooling. The objective function is the sum of the elements of the difference distance matrix between the two molecules generated by continual reordering of one molecule. Extensive tests of the algorithms have been performed on random coordinate data together with two related protein structures. Combinatorial problems, inherent in the assignment of atom correspondences, are effectively overcome by simulated annealing. The algorithms outlined here can readily optimize molecular matching problems with 150 atoms.     

Prediction of peptide conformation by multicanonical algorithm: New approach to the multiple-minima problem

We apply a recently developed method, the multicanonical algorithm, to the problem of tertiary structure prediction of peptides and proteins. As a simple example to test the effectiveness of the algorithm, metenkephalin is studied and the ergodicity problem, or multiple-minima problem, is shown to be overcome by this algorithm. The lowest-energy conformation obtained agrees with that determined by other efficient methods such as Monte Carlo simulated annealing. The superiority of the present method to simulated annealing lies in the fact that the relationship to the canonical ensemble remains exactly controlled. Once the multicanonical parameters are determined, only one simulation run is necessary to obtain the lowest-energy conformation and further the results of this one run can be used to calculate various thermodynamic quantities at any temperature. The latter point is demonstrated by the calculation of the average potential energy and specific heat as functions of temperature. © John Wiley & Sons, Inc.     

Simulated annealing to locate various stationary points in semiempirical methods

In this work, an algorithm was developed to study the potential energy surfaces in the coordinate spaces of molecules by a nonlocal way, in contrast to classic energy minimizers as the BFGS or the DFP method. This algorithm, based on the specificities of semiempirical methods, mixes simulated annealing and local searches to reduce computation costs. By this technique, the global energy minimum can be localized. Moreover, local minima that are close in energy to the global minimum are also obtained. If the search is not only for minima but for all stationary points (minima, saddle points…), then the energy is replaced by the gradient norm, which reaches its minimum values at stationary points. The annealing process is stopped before having accurately reached the global minimum and generates a list of geometries whose energies (respectively, whose gradients) are optimized by local minimizers. This list of geometries is shortened from the nearly equivalent geometries by a dynamic single-clustering analysis. The energy/gradient local minimizers act on the clustered list to produce a set of minima/stationary points. A targeted search of these points and reduction of the costs are reached by the way of several penalty functions. They eliminate—without energy calculation—most of the points generated by random walks on the potential energy surface. These penalty functions (on the total moment of inertia or on interatomic distances) are specific to the class of problem studied. They account for the nonrupture of either specific chemical bonds or rings in cyclic molecules, they assure that molecular systems are kept bonded, and they avoid the collapsing of atoms. © 1992 John Wiley & Sons, Inc.     

Genetic algorithm for the determination of binodal curves in ternary systems polymer-liquid(1)-liquid(2) and polymer(1)-polymer(2)-solvent

A simple genetic algorithm for the numerical evaluation of binodal curves in ternary systems polymer-liquid (1)-liquid (2) and polymer (1)-polymer (2)-solvent is presented. The technique exploits a specifically developed restarting technique based on a combined elitist and zooming strategy on the last population at each iteration. The objective function (fitness) is represented by the weighted sum of the squared differences of chemical potentials of the two phases of each component, obtained evaluating first derivatives of Gibbs free energy of the mixture with respect to the number of moles of the components. The method proposed (a) is numerically stable since it does not require the evaluation of first derivatives of the objective function and (b) can be applied in a wide range of cases changing the equation of state. Several comparisons with simplified iterative procedures presented in the past in the technical literature both for mixtures of two polymers with identical characteristics in a solvent and for mixtures of solvent-nonsolvent-polymer with solvent-polymer interaction parameter equal to zero are reported. Finally, a comparison between present results and the "alternating tangent approach" is reported for two technically meaningful binary systems, when a simplified PC-SAFT equation of state is adopted. The comparisons show that reliable results can be obtained by means of the algorithm proposed and suggest that the procedure presented can be used for practical purposes.     

Applications of simulated annealing to the conformational analysis of flexible molecules

We describe in this article our solution to the global minimum problem which uses the simulated annealing algorithm of Kirkpatrick. This method is a Metropolis (e-ΔE/kT) Monte Carlo sampling of conformation space with simultaneous constraint of the search by lowering the temperature T so that the search converges on the global minimum. The Anneal-Conformer program has been extensively tested with peptides and organic molecules using either the Amber or MM2 force fields. A history file of the simulated annealing process allows reconstruction of the random walk in conformation space for subsequent examination. Thus plots of distance and dihedral angle changes during the search for the global minimum can be examined to deduce molecular shape and flexibility. A separate program Conf-Gen reads the history file and extracts all low energy conformations visited during the run.     

Cluster algorithm to perform parallel Monte Carlo simulation of atomistic systems

We propose an efficient algorithm to perform Monte Carlo simulations of dense systems using multiple particle moves. The method is intended to be used in the atomistic simulation of complex systems, where the computing requirements for a single simulation run make advisable the use of parallel computing. The algorithm is based on the use of steps in which all the particle positions of the system are perturbed simultaneously. A division of the system in clusters of particles is performed, using a bonding criterion which makes feasible that the acceptance or rejection of the new particle coordinates can be carried out independently for each cluster.     

Encoding of polycyclic Si-containing molecules for determining species uniqueness in automated mechanism generation

Automated mechanism generation is an attractive way to understand the fundamental kinetics of complex reaction systems such as silicon hydride clustering chemistry. It relies on being able to tell molecules apart as they are generated. The graph theoretic foundation allows molecules to be identified using unique notations created from their connectivity. To apply this technique to silicon hydride clustering chemistry, a molecule canonicalization and encoding algorithm was developed to handle complex polycyclic, nonplanar species. The algorithm combines the concepts of extended connectivity and the idea of breaking ties to encode highly symmetric molecules. The connected components in the molecules are encoded separately and reassembled using a depth-first search method to obtain the correct string codes. A revised cycle-finding algorithm was also developed to properly select the cycles used for ring corrections when thermodynamic properties were calculated using group additivity. In this algorithm, the molecules are expressed explicitly as trees, and all linearly independent cycles of every size in the molecule are found. The cycles are then sorted according to their size and functionality, and the cycles with higher priorities will be used to include ring corrections. Applying this algorithm, more appropriate cycle selection and more accurate estimation of thermochemical properties of the molecules can be obtained.     

A reliable prediction of the global phase stability for liquid–liquid equilibrium through the simulated annealing algorithm: Application to NRTL and UNIQUAC equations

The simulated annealing algorithm is introduced to search the global optimal solutions for the multipeak phenomena which generally exist in the phase stability problems with continuous variables. The Gibbs free energy criterion was modeled by the NRTL and UNIQUAC activity coefficient equations. When previous approaches fail, it is usually because they locate local minima due to the nonconvex and nonlinear natures of the models used to predict phase equilibrium. In this paper, the preliminary results show that the global minimum of the tangent plane distance function (TPDF) can be obtained by using the simulated annealing algorithm. The effects of the initial and the final values of the control parameter, the decrement of the control parameter and the length of the Markov chains are analyzed. The optimal `cooling schedule' was obtained according to the calculation results of the phase stability problems for one ternary mixture. The liquid–liquid equilibrium compositions were calculated by the Newton–Raphson method on the basis of the global minimum of TPDF. The results of four examples show that the simulated annealing algorithm can effectively solve the global phase stability problem.     

Molecular structure matching by simulated annealing. II. An exploration of the evolution of configuration landscape problems

Summary This paper considers some of the landscape problems encountered in matching molecules by simulated annealing. Although the method is in theory ergodic, the global minimum in the objective function is not always encountered. Factors inherent in the molecular data that lead the trajectory of the minimization away from its optimal route are analysed. Segments comprised of the C atoms of dihydrofolate reductase are used as test data. The evolution of a reverse ordering landscape problem is examined in detail. Where such patterns in the data could lead to incorrect matches, the problem can in part be circumvented by assigning an initial random ordering to the molecules.     

A new hybrid algorithm for finding the lowest minima of potential surfaces: approach and application to peptides

A new algorithm is presented for finding the global minimum, and other low-lying minima, of a potential energy surface (PES) of biological molecules. The algorithm synergetically combines three well-known global optimization methods: the diffusion equation method (DEM), which involves smoothing the PES; a simulated annealing (SA) algorithm; and evolutionary programming (EP), whose population-oriented approach allows for a parallel search over different regions of the PES. Tests on five peptides having between 6 and 9 residues show that the code implementing the new combined algorithm is efficient and is found to outperform the constituent methods, DEM and SA. Results of the algorithm, in the gas phase and with the GBSA implicit solvent model, are compared with crystallographic data for the test peptides; good accord is found in all cases. Also, for all but one of the examples, our hybrid algorithm finds a minimum deeper than those obtained by a very extensive scan. TINKERs implementation of the OPLS-AA force field is employed for the structure prediction. The results show that the new algorithm is a powerful structure predictor, when a reliable potential function is available. Our implementation of the algorithm is time-efficient, and requires only modest computational resources. Work is underway on applications of the new algorithm to structural prediction of proteins and other biological macro-molecules.     

Development and validation of an improved algorithm for overlaying flexible molecules

A program for overlaying multiple flexible molecules has been developed. Candidate overlays are generated by a novel fingerprint algorithm, scored on three objective functions (union volume, hydrogen-bond match, and hydrophobic match), and ranked by constrained Pareto ranking. A diverse subset of the best ranked solutions is chosen using an overlay-dissimilarity metric. If necessary, the solutions can be optimised. A multi-objective genetic algorithm can be used to find additional overlays with a given mapping of chemical features but different ligand conformations. The fingerprint algorithm may also be used to produce constrained overlays, in which user-specified chemical groups are forced to be superimposed. The program has been tested on several sets of ligands, for each of which the true overlay is known from protein–ligand crystal structures. Both objective and subjective success criteria indicate that good results are obtained on the majority of these sets.     

Monte Carlo free energy calculations using electronic structure methods

The molecular mechanics-based importance sampling function (MMBIF) algorithm [R. Iftimie, D. Salahub, D. Wei, and J. Schofield, J. Chem. Phys. 113, 4852 (2000)] is extended to incorporate semiempirical electronic structure methods in the secondary Markov chain, creating a fully quantum mechanical Monte Carlo sampling method for simulations of reactive chemical systems which, unlike the MMBIF algorithm, does not require the generation of a system-specific force field. The algorithm is applied to calculating the potential of mean force for the isomerization reaction of HCN using thermodynamic integration. Constraints are implemented in the sampling using a modification of the SHAKE algorithm, including that of a fixed, arbitrary reaction coordinate. Simulation results show that sampling efficiency with the semiempirical secondary potential is often comparable in quality to force fields constructed using the methods suggested in the original MMBIF work. The semiempirical based importance sampling method presented here is a useful alternative to MMBIF sampling as it can be applied to systems for which no suitable MM force field can be constructed.