Restricted Hartree–Fock approximation. II. Computational aspects of the direct minimization procedure

A RHF energy minimization procedure based on the treatment outlined in Part I of this series of articles is presented. Test calculations performed on several closed- and open-shell systems show that the present procedure is definitely superior to the conventional SCF methods. In particular, the convergence of this procedure is ensured, the rate of convergency is high, and the computational cost of each cycle is low.     

Validated computing approach for high-pressure chemical and multiphase equilibrium

For the computation of chemical and phase equilibrium at constant temperature and pressure, there have been proposed a wide variety of problem formulations and numerical solution procedures, involving both direct minimization of the Gibbs energy and the solution of equivalent nonlinear equation systems. Still, with very few exceptions, these methodologies may fail to solve the chemical and phase equilibrium problem correctly. Nevertheless, there are many existing solution methods that are extremely reliable in general and fail only occasionally. To take good advantage of this wealth of available techniques, we demonstrate here an approach in which such techniques can be combined with procedures that have the power to validate results that are correct, and to identify results that are incorrect. Furthermore, in the latter case, corrective feedback can be provided until a result that can be validated as correct is found. The validation procedure is deterministic, and provides a mathematical and computational guarantee that the global minimum in the Gibbs energy has been found. To demonstrate this validated computing approach to the chemical and phase equilibrium problem, we present several examples involving reactive and nonreactive components at high pressure, using cubic equation-of-state models.     

Protein structure prediction using a combination of sequence homology and global energy minimization I. Global energy minimization of surface loops

A procedure has been developed for global energy minimization of surface loops of proteins in the presence of a fixed core. The ECEPP potential function has been modified to allow more accurate representations of hydrogen bond interactions and intrinsic torsional energies. A computationally efficient representation of hydration free energy has been introduced. A local minimization procedure has been developed that uses a cutoff distance, minimization with respect to subsets of degrees of freedom, analytical second derivatives, and distance constraints between rigid segments to achieve efficiency in applications to surface loops. Efficient procedures have been developed for deforming segments of the initial backbone structure and for removing overlaps. Global energy minimization of a surface loop is accomplished by generating a sequence (or a trajectory) of local minima, the component steps of which are generated by searching collections of local minima obtained by deforming seven-residue segments of the surface loop. The search at each component step consists of the following calculations: (1) A large collection of backbone structures is generated by deforming a seven-residue segment of the initial backbone structure. (2) A collection of low-energy backbone structures is generated by applying local energy minimization to the resulting collection of backbone structures (interactions involving side chains that will be searched in this component step are not included in the energy). (3) One low-energy side-chain structure is generated for each of the resulting low-energy backbone structures. (4) A collection of low-energy local minima is generated by applying local energy minimization to the resulting collection of structures. (5) The local minimum with the lowest energy is retained as the next point of the trajectory. Applications of our global search procedure to surface segments of bovine pancreatic trypsin inhibitor (BPTI) and bovine trypsin suggest that component-step searches are reasonably complete. The computational efficiency of component-step searches is such that trajectories consisting of about 10 component steps are feasible using an FPS-5200 array processor. Our procedure for global energy minimization of surface loops is being used to identify and correct problems with the potential function and to calculate protein structure using a combination of sequence homology and global energy minimization.     

Pattern recognition in the prediction of protein structure. II. Chain conformation from a probability-directed search procedure

A procedure that finds the most probable conformational states of a protein chain is described. Single-residue conformations are represented in terms of four conformational states, α, ?, α*, and ?*. The conformation of the entire chain is represented by a sequence of single-residue conformational states; the distinct conformations in this representation are called “chain-states.” The first article in this series described a procedure that computes tripeptide conformational probabilities from the amino acid sequence using pattern recognition techniques. The procedure described in this article uses the tripeptide probabilities to estimate the probabilities of the chain-states. The chain-state probability estimator is a product of conditional and marginal probabilities (obtained from the tripeptide probabilities), with a penalty factor to eliminate conformations containing α-helices and ?-strands of excessive length. The probability estimator considers short-range conformational information, medium-range sequence information and some simple long-range information (through the restrictions on helix and strand lengths). Energy minimization calculations can be carried out in the region of conformational space corresponding to a particular chain-state. By selecting the most probable chain-states, the search can be focused on the most probable, or “important,” regions of the conformational space. These energy calculations are described in the third article of the series. The complete procedure described by the three articles is called PRISM, for pattern recognition-based importance sampling minimization.     

Developing the molecular modelling of diffusion in zeolites as a high throughput catalyst screening technique

Molecular modelling techniques have been used to screen zeolite catalysts for their suitability for organic synthesis. For example, we have used these techniques to study the alkylation of aromatic molecules which are important in the fine-chemical and drug industries. A survey of all such efforts is reviewed in this article. The application of molecular modelling techniques in a systematic manner is an efficient first step in the design of zeolite catalysts. As a qualitative screening tool, molecular graphics is used to visualize how well the reactant and product molecules fit inside the pores of the zeolites. Using a hybrid of several molecular modelling methods, which combines molecular dynamics (MD) and Monte Carlo methods with energy minimization, it is possible to determine the minimum energy locations of the molecules inside the zeolites cages. The minimum energy configurations determined by this hybrid method are taken as a starting point for diffusion of the molecules through the zeolite channels. When a molecule is allowed to diffuse through zeolite channel, the molecule attains some maxima and minima in its diffusion energy profile. From the differences between a maximum and a minimum energy configuration, the diffusion energy barrier for the molecule can be calculated in the zeolites. By comparing the diffusion energy barriers for various isomers of a molecule in different zeolites, it is possible to find out the most suitable zeolite for achieving the required shape-selectivity. In addition, factors influencing the diffusivity of the molecules and consequently the shape selectivity are derived. The list of factors and their relative importance are analysed to derive valuable guidelines to design shape-selective zeolite catalysts for a given reaction. Thus, the ultimate aim of these studies is to develop a high throughput computational screening process for the selection of shape-selective zeolite catalysts for various reactions. The dynamic behaviour of molecules inside the pores of zeolites can be studied using MD methods. Since MD is computationally time consuming, it is more efficient to screen the possible zeolite catalysts by energy minimization methods and then perform MD in specific zeolites. More accurate values of diffusivity of the molecules can be calculated using MD methods, and these values can be correlated with the shape-selectivity observed experimentally and /or derived from diffusion energy barrier calculations.     

Pattern recognition in the prediction of protein structure. III. An importance-sampling minimization procedure

A procedure that generates random conformations of a protein chain, and then applies energy minimization to find the structure of lowest energy, is described. Single-residue conformations are represented in terms of four conformational states, α, ?, α*, and ?*. Each state corresponds to a rectangular region in the ?, ψ map. The conformation of an entire chain is then represented by a sequence of single-residue conformational states. The distinct “chain-states” in this representation correspond to multidimensional rectangular regions in the conformational space of the whole protein. A set of highly-probable chain-states can be predicted from the amino acid sequence using the pattern recognition procedure developed in the first two articles of this series. The importance-sampling minimization procedure of the present article is then used to explore the regions of conformational space corresponding to each of these chain-states. The importance-sampling procedure generates a number of random conformations within a particular multidimensional rectangular region, sampling most densely from the most probable, or “important,” sections of the ?, ψ map. All values of ? and ψ are allowed, but the less-probable values are sampled less often. To achieve this, the random values of ? and Φ are generated from bivariate gaussian distributions that are determined from known X-ray structures. Separate gaussian distributions are used for proline residues in the α and ? states, for glycine residues in the α, ?, α*, and ?* states, and for ordinary residues involved in 29 different tripeptide conformations. Energy minimization is then applied to the randomly-generated structures to optimize interactions and to improve packing. The final energy values are used to select the best structures. The importance-sampling minimization procedure is tested on the avian pancreatic polypeptide, using chain-states predicted from the amino acid sequence. The conformation having the lowest energy is very similar to the X-ray conformation.     

Gradient gravitational search: An efficient metaheuristic algorithm for global optimization

The adaptation of novel techniques developed in the field of computational chemistry to solve the concerned problems for large and flexible molecules is taking the center stage with regard to efficient algorithm, computational cost and accuracy. In this article, the gradient‐based gravitational search (GGS) algorithm, using analytical gradients for a fast minimization to the next local minimum has been reported. Its efficiency as metaheuristic approach has also been compared with Gradient Tabu Search and others like: Gravitational Search, Cuckoo Search, and Back Tracking Search algorithms for global optimization. Moreover, the GGS approach has also been applied to computational chemistry problems for finding the minimal value potential energy of two‐dimensional and three‐dimensional off‐lattice protein models. The simulation results reveal the relative stability and physical accuracy of protein models with efficient computational cost. © 2015 Wiley Periodicals, Inc.     

Efficiency of monte carlo minimization procedures and their use in analysis of NMR data obtained from flexible peptides

The Monte Carlo minimization (MCM) method of Li and Scheraga is an efficient tool for generating low energy minimized structures of peptides, in particular the global energy minimum (GEM). In a recent article we proposed an enhancement to MCM, called the free energy Monte Carlo minimization (FMCM) procedure. With FMCM the conformational search is carried out with respect to the harmonic free energy, which approximates the free energy of the potential energy wells around the energy minimized structures (these wells are called localized microstates). In this work we apply both methods to the pentapeptide Leu-enkephalin described by the potential energy function ECEPP, and study their efficiency in identifying the GEM structure as well as the global harmonic free energy (GFM) structure. We also investigate the efficiency of these methods to generate localized microstates, which pertain to different energy and harmonic free energy intervals above the GEM and GFM, respectively. Such microstates constitute an important ingredient of our statistical mechanical methodology for analyzing nuclear magnetic resonance data of flexible peptides. Aspects of this methodology related to the stability properties of the localized microstates are examined. © 1997 by John Wiley & Sons, Inc.     

Docking to RNA via root-mean-square-deviation-driven energy minimization with flexible ligands and flexible targets

Structure-based drug design is now well-established for proteins as a key first step in the lengthy process of developing new drugs. In many ways, RNA may be a better target to treat disease than a protein because it is upstream in the translation pathway, so inhibiting a single mRNA molecule could prevent the production of thousands of protein gene products. Virtual screening is often the starting point for structure-based drug design. However, computational docking of a small molecule to RNA seems to be more challenging than that to protein due to the higher intrinsic flexibility and highly charged structure of RNA. Previous attempts at docking to RNA showed the need for a new approach. We present here a novel algorithm using molecular simulation techniques to account for both nucleic acid and ligand flexibility. In this approach, with both the ligand and the receptor permitted some flexibility, they can bind one another via an induced fit, as the flexible ligand probes the surface of the receptor. A possible ligand can explore a low-energy path at the surface of the receptor by carrying out energy minimization with root-mean-square-distance constraints. Our procedure was tested on 57 RNA complexes (33 crystal and 24 NMR structures); this is the largest data set to date to reproduce experimental RNA binding poses. With our procedure, the lowest-energy conformations reproduced the experimental binding poses within an atomic root-mean-square deviation of 2.5 A for 74% of tested complexes.     

The effect of numerical error on the reproducibility of molecular geometry optimizations

Geometry optimization is one of the most often applied techniques in computational drug discovery. Although geometry optimization routines are generally deterministic, the minimization trajectories can be extremely sensitive to initial conditions, especially in case of larger systems such as proteins. Simple manipulations such as coordinate transformations (translations and rotations), file saving and retrieving, and hydrogen addition can introduce small variations ( approximately 0.001 A) in the starting coordinates which can drastically affect the minimization trajectory. With large systems, optimized geometry differences of up to 1 A RMSD and final energy differences of several kcal/mol can be observed when using many commercially available software packages. Differences in computer platforms can also lead to differences in minimization trajectories. Here we demonstrate how routine structure manipulations can introduce small variations in atomic coordinates, which upon geometry optimization, can give rise to unexpectedly large differences in optimized geometries and final energies. We also show how the same minimizations run on different computer platforms can also lead to different results. The implications of these findings on routine computational chemistry procedures are discussed.     

Ab initio calculations of optical absorption spectra: solution of the Bethe-Salpeter equation within density matrix perturbation theory

We describe an ab initio approach to compute the optical absorption spectra of molecules and solids, which is suitable for the study of large systems and gives access to spectra within a wide energy range. In this approach, the quantum Liouville equation is solved iteratively within first order perturbation theory, with a Hamiltonian containing a static self-energy operator. This procedure is equivalent to solving the statically screened Bethe-Salpeter equation. Explicit calculations of single particle excited states and inversion of dielectric matrices are avoided using techniques based on density functional perturbation theory. In this way, full absorption spectra may be obtained with a computational workload comparable to ground state Hartree-Fock calculations. We present results for small molecules, for the spectra of a 1 nm Si cluster in a wide energy range (20 eV), and for a dipeptide exhibiting charge transfer excitations.     

Quantum embedding for material chemistry based on domain separation and open subsystems

This perspective considers two theories we recently proposed to perform quantum embedding calculations for chemical systems: domain-separated density functional theory (DS-DFT) and locally coupled open subsystems (LCOS). The development includes both the fundamentals of each theory as well as potential applications, some technical aspects, and related challenges. DS-DFT is suited to study intramolecular effects, where one can apply a high level of theory (based on DFT or wave function theory) to a region of interest inside a molecule or solid and lower level theory elsewhere, with smooth switching between the regions. LCOS, in contrast, is a fragment-based embedding, which offers computational advantages to study intermolecular behavior such as electron hopping, spin-environment interaction, and charge-transfer excitations. However, both theories can exchange roles when appropriate. In addition, these theories allow for control of computational scaling of their algorithms. We explore paths to determine the charge-transfer operator used in LCOS, and suggest an auxiliary energy minimization that can provide a practical estimate to this operator. We also briefly discuss how to implement density fitting techniques in domain separation, and how domain separation can be used for pure wave function-based embedding.     

A powerful truncated Newton method for potential energy minimization

With advances in computer architecture and software, Newton methods are becoming not only feasible for large-scale nonlinear optimization problems, but also reliable, fast and efficient. Truncated Newton methods, in particular, are emerging as a versatile subclass. In this article we present a truncated Newton algorithm specifically developed for potential energy minimization. The method is globally convergent with local quadratic convergence. Its key ingredients are: (1) approximation of the Newton direction far away from local minima, (2) solution of the Newton equation iteratively by the linear Conjugate Gradient method, and (3) preconditioning of the Newton equation by the analytic second-derivative components of the “local” chemical interactions: bond length, bond angle and torsional potentials. Relaxation of the required accuracy of the Newton search direction diverts the minimization search away from regions where the function is nonconvex and towards physically interesting regions. The preconditioning strategy significantly accelerates the iterative solution for the Newton search direction, and therefore reduces the computation time for each iteration. With algorithmic variations, the truncated Newton method can be formulated so that storage and computational requirements are comparable to those of the nonlinear Conjugate Gradient method. As the convergence rate of nonlinear Conjugate Gradient methods is linear and performance less predictable, the application of the truncated Newton code to potential energy functions is promising.     

An analysis of the structural and energetic properties of deoxyribose by potential energy methods

We discuss the three fundamental issues of a computational approach in structure prediction by potential energy minimization, and analyze them for the nucleic acid component deoxyribose. Predicting the conformation of deoxyribose is important not only because of the molecule's central conformational role in the nucleotide backbone, but also because energetic and geometric discrepancies from experimental data have exposed some underlying uncertainties in potential energy calculations. The three fundamental issues examined here are: (i) choice of coordinate system to represent the molecular conformation; (ii) construction of the potential energy function; and (iii) choice of the minimization technique. For our study, we use the following combination. First, the molecular conformation is represented in cartesian coordinate space with the full set of degrees of freedom. This provides an opportunity for comparison with the pseudorotation approximation. Second, the potential energy function is constructed so that all the interactions other than the nonbonded terms are represented by polynomials of the coordinate variables. Third, two powerful Newton methods that are globally and quadratically convergent are implemented: Gill and Murray's Modified Newton method and a Truncated Newton method, specifically developed for potential energy minimization. These strategies have produced the two experimentally-observed structures of deoxyribose with geometric data (bond angles and dihedral angles) in very good agreement with experiment. More generally, the application of these modeling and minimization techniques to potential energy investigations is promising. The use of cartesian variables and polynomial representation of bond length, bond angle and torsional potentials promotes efficient second-derivative computation and, hence, application of Newton methods. The truncated Newton, in particular, is ideally suited for potential energy minimization not only because the storage and computational requirements of Newton methods are made manageable, but also because it contains an important algorithmic adaptive feature: the minimization search is diverted from regions where the function is nonconvex and is directed quickly toward physically interesting regions.     

Derivatives of molecular surface area and volume: simple and exact analytical formulas

The computational effort of biomolecular simulations can be significantly reduced by means of implicit solvent models in which the energy generally contains a correction depending on the surface area and/or the volume of the molecule. In this article, we present simple derivation of exact, easy-to-use analytical formulas for these quantities and their derivatives with respect to atomic coordinates. In addition, we provide an efficient, linear-scaling algorithm for the construction of the power diagram required for practical implementation of these formulas. Our approach is implemented in a C++ header-only template library.     

Correlation between the Estrada index and π-electronic energies for benzenoid hydrocarbons with applications to boron nanotubes

In this article, we present an efficient computer-based computational technique to compute the energy and Estrada index of graphs. It is shown that our computational method is more efficient and bears less computational and algorithmic complexity. We use our method to show the main result of this article, which asserts that the Estrada index correlates with the π-electronic energies of lower benzenoid hydrocarbons with correlation coefficient 0.9993. This enhances the practical applicability of the Estrada index and warrants its further usage in quantitative structure activity relationships. We further apply our computational technique in computing the energy and Estrada index of two infinite families of boron triangular nanotubes. We perform simulation based on certain computer software packages to study the graph-theoretic behavior of the obtained results. Our results help to correlate certain physicochemical properties of underlying chemical structures of these nanotubes.     

Monte Carlo computation of ground‐state energy derivatives

Monte Carlo is a simple technique, which uses random numbers to compute ground‐state energies of small molecules (and quantum systems in general). The results always have a small statistical error, which poses a major obstacle when estimating properties defined as ground‐state‐energy derivatives (such as the molecule's geometry, its vibrational frequencies, polarizabilities, etc.). In this article, we present and demonstrate an approach that makes an accurate Monte–Carlo estimation of such derivatives possible. This is achieved by realizing that the simulation constitutes an autocorrelated stochastic process, whose proper analysis then enables us to estimate various energy derivatives as a combination of total correlation between readily computable quantities. The resulting procedure is a natural extension of the usual Monte Carlo algorithm for computing the ground‐state energy, with relatively small computational overhead. © 2007 Wiley Periodicals, Inc. Int J Quantum Chem, 2008     

Monte Carlo simulation of RRKM unimolecular decomposition in molecular beam experiments. I. Basic considerations and calculational procedure

A simulation method is described for the comparison of a molecular decomposition theory, based on the fundamental RRKM theory, with crossed molecular beam experiments. In the present formulation, the method is applied to the case with long-lived collision complexes surrounded by a centrifugal barrier. The procedure uses Monte Carlo techniques to simulate the formation and decomposition of the complexes, thus effectively solving the high-dimensional integrals resulting but seldom solved in the analytical treatments. Several computational problems, like singularities in the c.m. angular distributions, have been circumvented in this procedure. Angular momentum conservation, the simultaneous interaction of several decomposition channels and strict flux conservation are included in the procedure. Dynamical features can be introduced, as will be shown in forthcoming papers.