Chapter 4 - Self-Consistent Field using Direct Inversion in Iterative Subspace Method and Quasi-Newton Vectors

A novel doubly iterative self-consistent field (SCF) approach using a semiempirical model Hamiltonian (denoted as the SMH algorithm) is proposed to accelerate the Hartree-Fock (HF) and density functional theory (DFT) calculations. This method first constructs the Fock matrix exactly in each SCF macroiteration, followed by a few SCF microiterations, where the Fock matrix is incrementally updated using an inexpensive semiempirical approximation. This leads to an improved wave function in each SCF macroiteration with minimal additional cost, and therefore a reduced number of exact Fock builds is required for SCF convergence. The SMH algorithm can be combined with conventional SCF convergence techniques such as level shifting, direct inversion in the iterative subspace (DIIS), and energy-DIIS (EDIIS). When integrated with DIIS, SMH enhances the convergence of simple organic molecules by approximately 10% compared to plain DIIS, with speedups of up to 60% for the more challenging transition metal systems compared to EDIIS + DIIS. Our results show that SMH is a reliable SCF accelerator that seldom deteriorates convergence and is highly robust.

The direct inversion in the iterative subspace (DIIS) method is applied to several simple SCF wave functions in an effective Fock matrix formulation. The following cases are treated: high-spin-restricted open shell, open-shell singlet, and two-configuration wave functions. Open-shell singlet states are described by a three-determinant 2×2 CAS expansion which is equivalent to Davidson’s nonorthogonal SCF method in the case of the first open-shell singlet. Very sharp convergence is usually obtained in less than 20 cycles. The method is applicable to slowly convergent or even inherently divergent cases, and able to enforce convergence to excited states not the lowest of their symmetry. For these simple wave functions, the present first order method is asymptotically more efficient than second-order methods. Examples are presented for H2O, H2O2, C2H4, F2, several states of NO2, C2H5, formaldehyde, and ketene.

For electronic structure calculations on large systems, solving the self-consistent-field (SCF) equations is one of the key bottlenecks. Density matrix search methods provide an efficient linear scaling approach for circumventing the traditional O(N3) diagonalization procedure for solving the SCF equations. The conjugate gradient density matrix search (CG-DMS) method is a successful implementation of this approach. An alternative density matrix search method, QN–DMS, employs direct inversion in the iterative subspace using a quasi-Newton (QN) step to estimate the error vector. For linear polyglycine chains of 10–100 residues, the present approach scales linearly and is 30% faster than CG-DMS. For clusters of up to 300 water molecules, this method shows smoother and efficient convergence, and displays nearly linear scaling.

We extend the application of the direct inversion in the iterative subspace (DIIS) technique to the ridge method for finding transition states (TS). The latter is not a quasi-Newton-type algorithm, which is the only class of geometry optimization methods that has been combined with DIIS. With this new combination, we obtain a factor of two speedup due to DIIS, similar to the DIIS-related speedup achieved in other methods including quasi-Newton geometry optimization and self-consistent field iterations. We also demonstrate that DIIS is useful even in cases where optimization is started far from the quadratic region of the TS, provided that only one previous iteration is used in the DIIS expansion. We compare the performance of the new ridge-DIIS method to that of the TS algorithm utilized in GAUSSIAN 92. We find that the computational cost of the former is similar (when both methods converge) to that of the latter. The examples considered in the paper include a novel TS found for an isomerization of a cluster of six Na atoms. Locating such a TS poses a known problem for second-derivatives-based algorithms that fail on very flat potential energy surfaces. Thus, the gradient-based ridge-DIIS method is the only TS search method that is robust, does not need second derivatives and/or an initial guess for the TS geometry, and whose performance matches or exceeds that of a second-derivatives-based algorithm.

Various geometry optimization techniques are systematically investigated. The rational function (RF) and direct inversion in the iterative subspace (DIIS) methods are compared and optimized for the purpose of geometry optimization. Various step restriction and line search procedures are tested. The model Hessian recently proposed by Lindh et al. has been used in conjunction with different Hessian update procedures. Optimization for over 30 molecules have been performed in Z-matrix coordinates, local normal coordinates, and curvilinear natural internal coordinates, using the same approximations for the Hessian in all cases. The most effective and stable procedure for optimization of equilibrium structures was found to be the DIIS minimization in natural internal coordinates using the BFGS update of the model Hessian. Our method shows faster overall convergence than all previously published methods for the same test suite of molecules. © 1997 John Wiley & Sons, Inc. J Comput Chem 18: 1473–1483, 1997

A geometry optimization method using an energy-represented direct inversion in the iterative subspace algorithm, GEDIIS, is introduced and compared with another DIIS formulation (controlled GDIIS) and the quasi-Newton rational function optimization (RFO) method. A hybrid technique that uses different methods at various stages of convergence is presented. A set of test molecules is optimized using the hybrid, GEDIIS, controlled GDIIS, and RFO methods. The hybrid method presented in this paper results in smooth, well-behaved optimization processes. The optimization speed is the fastest among the methods considered.

The efficacy of a method based on the direct inversion in the iterative subspace (DIIS) in accelerating the approach to self consistency in the calculation of theMCSCF wavefunction using a novel algorithm developed earlier, is compared with that of a simple damping technique. Although the ‘damping’ turns out to be ineffective in the ‘quadratic region’, it accelerates remarkably in the rate of descent on the energy hypersurface in the early stages of the iterative process which leads to an impressive overall increase in the rate of approach to self consistency. TheDIIS based procedure turns out to be ineffective when coupled to the present method and is plagued by ill conditioning problems. Calculations are done to compute the equilibrium geometrical parameter, charge density on different atoms, and dipole moment of HNO molecule in the lowest1,3 nπ* states at theINDO/2-MCSCF level.

The nonperturbative calculation of molecular properties in magnetic fields requires the evaluation of integrals over complex-valued Gaussian-type London atomic orbitals (LAOs). With these orbitals, the calculation of four-centre electron-repulsion integrals (ERIs) is particularly demanding, because their permutational symmetry is lowered, and because complex algebra is required. We have implemented the resolution-of-the-identity (RI) approximation for LAOs in the TURBOMOLE program package. With respect to LAOs, employing the RI approximation is particularly beneficial, because the auxiliary basis set may always be chosen to be real-valued. As a consequence, the two-centre integrals in the RI approximation remain real-valued, and the three-centre integrals possess the same permutational symmetry as their real-valued counterparts. Compared to a direct calculation of four-centre ERIs over LAOs, using the RI approximation thus not only reduces the scaling of the integral evaluation, but also increases the efficiency by an additional factor of at least two. By using other well-established methods such as Cauchy–Schwarz screening, the difference-density approach, and Pulay's direct inversion in the iterative subspace (DIIS), the efficiency of nonperturbative calculations in magnetic fields can be increased even further.

Sparse eigenproblems are important for various applications in computer graphics. The spectrum and eigenfunctions of the Laplace–Beltrami operator, for example, are fundamental for methods in shape analysis and mesh processing. The Subspace Iteration Method is a robust solver for these problems. In practice, however, Lanczos schemes are often faster. In this article, we introduce the Hierarchical Subspace Iteration Method (HSIM) , a novel solver for sparse eigenproblems that operates on a hierarchy of nested vector spaces. The hierarchy is constructed such that on the coarsest space all eigenpairs can be computed with a dense eigensolver. HSIM uses these eigenpairs as initialization and iterates from coarse to fine over the hierarchy. On each level, subspace iterations, initialized with the solution from the previous level, are used to approximate the eigenpairs. This approach substantially reduces the number of iterations needed on the finest grid compared to the non-hierarchical Subspace Iteration Method. Our experiments show that HSIM can solve Laplace–Beltrami eigenproblems on meshes faster than state-of-the-art methods based on Lanczos iterations, preconditioned conjugate gradients, and subspace iterations.

The commutator direct inversion of the iterative subspace (commutator DIIS or C-DIIS) method developed by Pulay is an efficient and the most widely used scheme in quantum chemistry to accelerate the convergence of self-consistent field (SCF) iterations in Hartree-Fock theory and Kohn-Sham density functional theory. The C-DIIS method requires the explicit storage of the density matrix, the Fock matrix, and the commutator matrix. Hence, the method can only be used for systems with a relatively small basis set, such as the Gaussian basis set. We develop a new method that enables the C-DIIS method to be efficiently employed in electronic structure calculations with a large basis set such as planewaves for the first time. The key ingredient is the projection of both the density matrix and the commutator matrix to an auxiliary matrix called the gauge-fixing matrix. The resulting projected commutator-DIIS method (PC-DIIS) only operates on matrices of the same dimension as that consists of Kohn-Sham orbitals. The cost of the method is comparable to that of standard charge mixing schemes used in large basis set calculations. The PC-DIIS method is gauge-invariant, which guarantees that its performance is invariant with respect to any unitary transformation of the Kohn-Sham orbitals. We demonstrate that the PC-DIIS method can be viewed as an extension of an iterative eigensolver for nonlinear problems. We use the PC-DIIS method for accelerating Kohn-Sham density functional theory calculations with hybrid exchange-correlation functionals, and demonstrate its superior performance compared to the commonly used nested two-level SCF iteration procedure. Furthermore, we demonstrate that in the context of ab initio molecular dynamics (MD) simulation with hybrid functionals one can extrapolate the gauge-fixing matrix to achieve the goal of extrapolating the entire density matrix implicitly along the MD trajectory. Numerical results indicate that the new method significantly reduces the number of SCF iterations per MD step, compared to the commonly used strategy of extrapolating the electron density.

Comparison of Lanczos and subspace iterations for hyperspherical reaction path calculations

The self-consistent field (SCF) procedure is the standard technique for solving the Hartree-Fock and Kohn-Sham density functional theory calculations, while convergence is not theoretically guaranteed. Direct minimization methods, such as the augmented Lagrangian method (ALM) and second-order SCF (SOSCF), obtain the SCF solution by minimizing the Lagrangian with the gradient. In SOSCF, molecular orbitals are optimized by truncating the Taylor expansion of a unitary matrix represented in exponential form to ensure the orthonormality condition. This study proposes an alternative algorithm for direct-energy minimization to obtain an SCF solution using ALM Lagrangian by adopting sequential Givens rotations between occupied and virtual orbitals. The Givens rotation corresponds to unitary transformations that guarantee orthogonality and avoid variational collapse. Complex gradients for sequential Givens rotation were obtained by the error back-propagation method, which is based on the chain rule. Illustrative applications clarified the features of the present DGR methods by comparing with other SCF algorithms such as direct inversion in iterative subspace, SOSCF, and ALM.

In this paper, a new wideband direction-of-arrival (DOA) estimation algorithm, which is an iterative coherent subspace method based on two-dimensional discrete Fourier transform (2D-DFT), is proposed. Two stages of focusing are accomplished in the algorithm. Firstly, a modified 2D-DFT, which has accomplished a preliminary focusing process, is exploited on the received data to achieve initial estimates that are more accurate than the results of beamforming methods. More importantly, in the first stage, the real array manifold is projected to a new virtual array manifold that is more suitable for subsequent focusing. Secondly, based on the transformed data, an iterative process in which each step works on the estimates available from previous iterations is applied to the transformed data to estimate the DOA. Because of the implementation of 2D-DFT, the proposed algorithm can achieve a higher resolution and a better convergence performance than conventional methods, especially when the incoming sources have irregular distributions in the frequency domain. Performance analysis and simulation results show that the proposed algorithm works properly and effectively.

Three algorithms for efficient solution of optimal control problems for high-dimensional systems are presented. Each bypasses the intermediate (and, unnecessary) step of open-loop model reduction. Each also bypasses the solution of the full Riccati equation corresponding to the LQR problem, which is numerically intractable for large n. Motivation for this effort comes from the field of model-based flow control, where open-loop model reduction often fails to capture the dynamics of interest (governed by the Navier–Stokes equation). Our minimum control energy method is a simplified expression for the well-known minimum-energy stabilizing control feedback that depends only on the left eigenvectors corresponding to the unstable eigenvalues of the system matrix A. Our Adjoint of the Direct-Adjoint method is based on the repeated iterative computation of the adjoint of a forward problem, itself defined to be the direct-adjoint vector pair associated with the LQR problem. Our oppositely-shifted subspace iteration (OSSI, the main new result of the present paper) method is based on our new subspace iteration method for computing the Schur vectors corresponding, notably, to the \(m\ll n\) central eigenvalues (near the imaginary axis) of the Hamiltonian matrix related to the Riccati equation of interest. Prototype OSSI implementations are tested on a low-order control problem to illustrate its behavior.

The subspace iteration method – Revisited