The trust-region self-consistent field method: towards a black-box optimization in Hartree-Fock and Kohn-Sham theories

The trust-region self-consistent field (TRSCF) method is presented for optimizing the total energy E(SCF) of Hartree-Fock theory and Kohn-Sham density-functional theory. In the TRSCF method, both the Fock/Kohn-Sham matrix diagonalization step to obtain a new density matrix and the step to determine the optimal density matrix in the subspace of the density matrices of the preceding diagonalization steps have been improved. The improvements follow from the recognition that local models to E(SCF) may be introduced by carrying out a Taylor expansion of the energy about the current density matrix. At the point of expansion, the local models have the same gradient as E(SCF) but only an approximate Hessian. The local models are therefore valid only in a restricted region-the trust region-and steps can only be taken with confidence within this region. By restricting the steps of the TRSCF model to be inside the trust region, a monotonic and significant reduction of the total energy is ensured in each iteration of the TRSCF method. Examples are given where the TRSCF method converges monotonically and smoothly, but where the standard DIIS method diverges.     

Extremely localized molecular orbitals (ELMO): a non-orthogonal Hartree-Fock method

A new optimization method for extremely localized molecular orbitals (ELMO) is derived in a non-orthogonal formalism. The method is based on a quasi Newton-Raphson algorithm in which an approximate diagonal-blocked Hessian matrix is calculated through the Fock matrix. The Hessian matrix inverse is updated at each iteration by a variable metric updating procedure to account for the intrinsically small coupling between the orbitals. The updated orbitals are obtained with approximately n ² operations. No n ³ processes such as matrix diagonalization, matrix multiplication or orbital orthogonalization are employed. The use of localized orbitals allows for the creation of high-quality initial “guess” orbitals from optimized molecular orbitals of small systems and thus reduces the number of iterations to converge. The delocalization effects are included by a Jacobi correction (JC) which allows the accurate calculation of the total energy with a limited number of operations. This extension, referred to as ELMO(JC), is a variational method that reproduces the Hartree-Fock (HF) energy with an error of less than 2 kcal/mol for a reduced total cost compared to standard HF methods. The small number of variables, even for a very large system, and the limited number of operations potentially makes ELMO a method of choice to study large systems. Received: 30 December 1996 / Accepted: 5 June 1997     

Methods for parallel computation of SCF NMR chemical shifts by GIAO method: Efficient integral calculation,multi-Fock algorithm,and pseudodiagonalization

We implemented our gauge-including atomic orbital (GIAO) NMR chemical shielding program on a workstation cluster, using the parallel virtual machine (PVM) message-passing system. On a modest number of nodes, we achieved close to linear speedup. This program is characterized by several novel features. It uses the new integral program of Wolinski that calculates integrals in vectorized batches, increases efficiency, and simplifies parallelization. The self-consistent field (SCF) step includes a multi-Fock algorithm, i.e., the simultaneous calculation of several Fock matrices with the same integral set, increasing the efficiency of the direct SCF procedure. The SCF diagonalization step, which is difficult to parallelize, has been replaced by pseudodiagonalization. The latter, widely used in semiempirical programs, becomes important in ab initio type calculations above a certain size, because the ultimate scaling of the diagonalization step is steeper than that of integral computation. Examples of the calculation of the NMR shieldings in large systems at the SCF level are shown. Parallelization of the density functional code is underway. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18: 816–825, 1997     

Accurate and efficient determination of higher roots in diagonalization of large matrices based in function restricted optimization algorithms

The problem of large‐scale matrix diagonalization is analyzed in the context of normal function optimization techniques with particular emphasis on the problem of obtaining high roots. New methods based on function restricted optimization algorithms are presented. The efficiency of these methods is illustrated for lowest and higher and degenerate roots of selected matrices. The diagonalization process is commonly carried out in a subspace, and involves a sort of optimization process, and the dimension of this subspace increases at each iteration. In addition, the success of a diagonalization method in obtaining a desired root strongly depends on the particular optimization procedure chosen. In this work, a rational function optimization procedure is presented that permits obtaining the lowest and higher eigenpairs in an efficient way. Update Hessian matrices formulae, routinely used in normal function optimization problems, are explored in the framework of diagonalization techniques. Finally, a diagonalization method with a fixed subspace dimension during the iterative process is presented. Some examples focused in lowest, higher and degenerate eigenpairs are discussed. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1375–1386, 2000     

Development of parallel density functional program using distributed matrix to calculate all-electron canonical wavefunction of large molecules

We developed a new parallel density-functional canonical molecular-orbital program for large molecules based on the resolution of the identity method. In this study, all huge matrices were decomposed and saved to the distributed local memory. The routines of the analytical molecular integrals and numerical integrals of the exchange-correlation terms were parallelized using the single program multiple data method. A conventional linear algebra matrix library, ScaLAPACK, was used for matrix operations, such as diagonalization, multiplication, and inversion. Anderson's mixing method was adopted to accelerate the self-consistent field (SCF) convergence. Using this program, we calculated the canonical wavefunctions of a 306-residue protein, insulin hexamer (26,790 orbitals), and a 133-residue protein, interleukin (11,909 orbitals) by the direct-SCF method. In regard to insulin hexamer, the total parallelization efficiency of the first SCF iteration was estimated to be 82% using 64 Itanium 2 processors connected at 3.2 GB/s (SGI Altix3700), and the calculation successfully converged at the 17-th SCF iteration. By adopting the update method, the computational time of the first and the final SCF loops was 229 min and 156 min, respectively. The whole computational time including the calculation before the SCF loop was 2 days and 17 h. This study put the calculations of the canonical wavefunction of 30,000 orbitals to practical use.     

Another way to implement the Powell formula for updating Hessian matrices related to transition structures

A way to update the Hessian matrix according to the Powell formula is given. With this formula one does not need to store the full Hessian matrix at any iteration. A method to find transition structures, which is a combination of the quasiNewton–Raphson augmented Hessian algorithm with the proposed Powell update scheme, is also given. The diagonalization of the augmented Hessian matrix is carried out by Lanczoslike methods. In this way, during all the optimization process, one avoids to store full matrices.     

A reformulation of the coupled perturbed self-consistent field equations entirely within a local atomic orbital density matrix-based scheme

To exploit the exponential decay found in numerical studies for the density matrix and its derivative with respect to nuclear displacements, we reformulate the coupled perturbed self-consistent field (CPSCF) equations and a quadratically convergent SCF (QCSCF) method for Hartree-Fock and density functional theory within a local density matrix-based scheme. Our D-CPSCF (density matrix-based CPSCF) and D-QCSCF schemes open the way for exploiting sparsity and to achieve asymptotically linear scaling of computational complexity with molecular size (M), in case of D-CPSCF for all (M) derivative densities. Furthermore, these methods are even for small molecules strongly competitive to conventional algorithms.     

Spinor optimization for a relativistic spin-dependent CASSCF program

Detailed formulae for the implementation of the multi-configuration SCF spinor optimization in a basis of Kramers pair 2-spinors – i.e. exploiting time-reversal symmetry – are presented. Full expressions for the spinor gradient and spinor Hessian elements are given in abstract form as well as within the usual CASSCF subspace division. As far as possible, the resulting terms are grouped to relativistic inactive and active Fock matrices, which have been introduced previously. Approximations for the Hessian are introduced so as to initialize it in an inverse Hessian update algorithm for a diagonal first approximation within the standard quasi-Newton-Raphson procedure. The effects of double group symmetry arising from spin dependence on Fock matrices and therefore gradient and Hessian are discussed and a group scheme for the implementation is proposed. Received: 14 January 1997 / Accepted: 3 February 1997     

Implementation of the diagonalization‐free algorithm in the self‐consistent field procedure within the four‐component relativistic scheme

A recently developed Thouless‐expansion‐based diagonalization‐free approach for improving the efficiency of self‐consistent field (SCF) methods (Noga and ?imunek, J. Chem. Theory Comput. 2010, 6, 2706) has been adapted to the four‐component relativistic scheme and implemented within the program package ReSpect. In addition to the implementation, the method has been thoroughly analyzed, particularly with respect to cases for which it is difficult or computationally expensive to find a good initial guess. Based on this analysis, several modifications of the original algorithm, refining its stability and efficiency, are proposed. To demonstrate the robustness and efficiency of the improved algorithm, we present the results of four‐component diagonalization‐free SCF calculations on several heavy‐metal complexes, the largest of which contains more than 80 atoms (about 6000 4‐spinor basis functions). The diagonalization‐free procedure is about twice as fast as the corresponding diagonalization. © 2014 Wiley Periodicals, Inc.     

On the restricted step method coupled with the augmented Hessian for the search of stationary points of any continuous function

In any optimization using the augmented Hessian technique, the step is not restricted to any length. Since the restriction of the step at each iteration is very important in order to achieve good convergence, we present a coupled method such that the augmented Hessian automatically gives both the adequate length of the step and the correct Hessian structure. The method is showed for the minima and saddle points of any order. © 1997 John Wiley & Sons, Inc.     

An efficient self-consistent field method for large systems of weakly interacting components

An efficient method for removing the self-consistent field (SCF) diagonalization bottleneck is proposed for systems of weakly interacting components. The method is based on the equations of the locally projected SCF for molecular interactions (SCF MI) which utilize absolutely localized nonorthogonal molecular orbitals expanded in local subsets of the atomic basis set. A generalization of direct inversion in the iterative subspace for nonorthogonal molecular orbitals is formulated to increase the rate of convergence of the SCF MI equations. Single Roothaan step perturbative corrections are developed to improve the accuracy of the SCF MI energies. The resulting energies closely reproduce the conventional SCF energy. Extensive test calculations are performed on water clusters up to several hundred molecules. Compared to conventional SCF, speedups of the order of (N/O)2 have been achieved for the diagonalization step, where N is the size of the atomic orbital basis, and O is the number of occupied molecular orbitals.     

Extended hessian matrix along the reaction coordinate

We shall prove by the coordinate-free covariant treatment that the vibrational frequency along the reaction coordinate is exactly given by the extended Hessian matrix in a rigged configuration space where the reaction coordinate is identified with a straight line. The diagonalization of the extended Hessian matrix is straightforward with no recourse to the prescribed procedure of eliminating the coupling matrix elements between the promoting mode of chemical reaction and the vibrational modes.     

A new direct CI method for large CI expansions in a small orbital space

A new direct CI method is presented, which is particularly suited for large CI expansions in a small orbital space. These are the type of expansions which are common in the CAS SCF method. Only one-electron coupling coefficients are stored, which leads to reduced elapsed times and storage requirements compared to earlier approaches. The two-electron coupling coefficients are implicitly created in the diagonalization step. The algorithm for updating the CI vector is formulated as the trace of a product of three matrices, Tr(A · D · I). By ordering the one-electron coupling coefficients (A) in a certain way the matrix D is easilly created as a sparse scalar product between these coefficients and the trial CI vector. The main computational step is then a simple matrix multiplication between the matrix D and the symmetry blocked integral matrix (1). This operation vectorizes very well on most vector processors. Another sparse scalar product between the resultant matrix and the coupling coefficients leads to the update of the CI coefficients. In a calculation on CRAY-1 with 30700 configurations, the two-electron part in a CI iteration required 10 s of which half went into the handling of the one-electron formula tape.     

Extracting effective normal modes from equilibrium dynamics at finite temperature

A general method for obtaining effective normal modes of a molecular system from molecular dynamics simulations is presented. The method is based on a localization criterion for the Fourier transformed velocity time-correlation functions of the effective modes. For a given choice of the localization function used, the method becomes equivalent to the principal mode analysis (PMA) based on covariance matrix diagonalization. On the other hand, a proper choice of the localization function leads to a novel method with a strong analogy with the usual normal mode analysis of equilibrium structures, where the Hessian system at the minimum energy structure is replaced by the thermal averaged Hessian, although the Hessian itself is never actually calculated. This method does not introduce any extra numerical cost during the simulation and bears the same simplicity as PMA itself. It can thus be readily applied to ab initio molecular dynamics simulations. Three such examples are provided here. First we recover effective normal modes of an isolated formaldehyde molecule computed at 20 K in very good agreement with the results of a normal mode analysis performed at its equilibrium structure. We then illustrate the applicability of the method for liquid phase studies. The effective normal modes of a water molecule in liquid water and of a uracil molecule in aqueous solution can be extracted from ab initio molecular dynamics simulations of these two systems at 300 K.     

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.     

Normal Coordinate Analysis for Polymer Systems: Capabilities and New Opportunities

Normal coordinate analysis is an important tool in studying the structure, dynamics, and physical properties of polymer systems. In this article the capabilities of normal coordinate analysis (NCA) are explored in some detail. The use of the eigenvalues and eigenvectors from NCA is catalogued for a wide variety of purposes: for assigning or interpreting polymer spectra, for structural determination, for constructing force fields, for computing heat capacity and other thermodynamic properties, and for computing other physical properties. Examples are given for crystals, melts, and amorphous systems. Also described are methods for characterizing the normal mode vectors that are especially useful for larger systems, in which a large amount of data must be analyzed or where visualization or animation fails. Finally, a recently developed method for eliminating negative eigenvalues in systems with tens of thousands of atoms, trajectory averaging, is presented. Also described are several advances in numerical linear algebra for speeding up the diagonalization phase and for computing physical properties without requiring full diagonalization of the Hessian matrix.     

An exact reformulation of the diagonalization step in electronic structure calculations as a set of second order nonlinear equations

A new formulation of the diagonalization step in self-consistent-field (SCF) electronic structure calculations is presented. It exactly replaces the diagonalization of the effective Hamiltonian with the solution of a set of second order nonlinear equations. The density matrix and/or the new set of occupied orbitals can be directly obtained from the resulting solution. This formulation may offer interesting possibilities for new approaches to efficient SCF calculations. The working equations can be derived either from energy minimization with respect to a Cayley-type parametrization of a unitary matrix, or from a similarity transformation approach.     

Integrating steepest-descent reaction pathways for large molecules

Exploring potential energy surfaces of large molecular systems can be quite challenging due to the increased number of nuclear degrees of freedom. Many techniques that are well-suited for small and moderate size systems require diagonalization of the energy second-derivative matrix. Since the cost of this step scales as O(N(atoms)(3)) (where N(atoms) is the number of atomic centers), such methods quickly become infeasible and are eventually rendered cost prohibitive. In this work, the recently developed Euler-based predictor-corrector reaction path integration method [H. P. Hratchian, M. J. Frisch, and H. B. Schlegel, J. Chem. Phys. 133, 224101 (2010)] is enhanced and proposed as a useful alternative to conventional reaction path following schemes in studies on very large systems. Because this integrator does not require Hessian diagonalization, the O(N(atoms)(3)) bottleneck afflicting other approaches is completely avoided. The effectiveness of the integrator in large system studies is demonstrated with an enzyme-catalyzed reaction employing an ONIOM (QM:MM) model chemistry and involving 5368 atomic centers.     

How does it become possible to treat delocalized and/or open-shell systems in fragmentation-based linear-scaling electronic structure calculations? The case of the divide-and-conquer method

The authors have developed a fragmentation-based linear-scaling electronic structure calculation strategy named the divide-and-conquer (DC) method, which has been implemented into the Gamess program package. Although there are many sorts of fragmentation-based linear-scaling schemes, most of them require the charge and spin multiplicity of each fragment a priori. Therefore, their applications to delocalized and/or open-shell systems have been limited. However, the DC method is a notable exception because the distribution of electrons in the entire system is automatically determined by the universal Fermi level. In this perspective, the authors have summarized the performance of the linear-scaling self-consistent field (SCF) and post-SCF calculations of delocalized and/or open-shell systems based on the DC method. Furthermore, some future prospects of the method have been discussed.     

Use of exchange maximization to generate starting vectors for self-consistent field calculations on metal cluster/adsorbate systems

Localized molecular orbitals (LMOs) derived from exchange maximization with respect to all atom-centered basis functions in the basis set are shown to generate a good starting electronic field for self-consistent field calculations on extended systems such as metal clusters, for which well-defined chemical bonds are not present. Examples studied are a cluster of 20 Ni atoms and the Pt(97)CO, Ag(43)/H(3)CNON, Ag(91)/H(2)CO, and vinylidene/Ni metal cluster plus adsorbate systems. It is also shown that improved starting vectors can be obtained by remixing a subset of the LMOs with the largest exchange eigenvalues through diagonalization of the Fock matrix computed with a null electronic field. Employing only a subset of the exchange-maximized LMOs in the first iterations, and then gradually expanding the space in which the diagonalizations are carried out in succeeding cycles, is shown to be an effective means of guiding the SCF procedure to the converged full-basis solution.