Direct Givens rotation method based on error back-propagation algorithm for self-consistent field solution.

Because the Hartree-Fock (HF) equations are non-linear in the molecular orbitals (MO’s), they must be solved iteratively to obtain a self-consistent set of MO’s (i.e. the SCF method). When iterations of the SCF process converge, the energy is stationary with respect to infinitesmal variations in the orbitals. Whether this stationary point is stable (a local minimum) or unstable (a maximum in one or more directions) depends on the second derivatives of the energy with respect to the orbital variations. The internal and external stabilities or instabilities of various constrained HF wavefunctions (real and complex, restricted, unrestricted and general HF) are discussed. A simple iterative scheme does not always lead to convergence of the SCF equations. The connection between SCF stability and convergence is examined and a number of methods for improving SCF convergence are discussed, including damping, level shifting, extrapolation, univariate search, direct inversion of iterative subspace (DIIS) and quadratically convergent SCF. Techniques for constructing better initial guess MO’s are also considered.

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.

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.

KSSOLV (Kohn-Sham Solver) is a MATLAB (Matrix Laboratory) tool-box for solving the Kohn-Sham density functional theory (KS-DFT) with the plane-wave basis set. In the KS-DFT calculations, the most expensive part is commonly the diagonalization of Kohn-Sham Hamiltonian in the self-consistent field (SCF) scheme. To enable a personal computer to perform medium-sized KS-DFT calculations that contain hundreds of atoms, we present a hybrid CPU-GPU implementation to accelerate the iterative diagonalization algorithms implemented in KSSOLV by using the MATLAB built-in Parallel Computing Toolbox. We compare the performance of KSSOLV-GPU on three types of GPU, including RTX3090, V100, and A100, with conventional CPU implementation of KSSOLV respectively and numerical results demonstrate that hybrid CPU-GPU implementation can achieve a speedup of about 10 times compared with sequential CPU calculations for bulk silicon systems containing up to 128 atoms.

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 present a subspace projection technique to conduct large-scale Kohn-Sham density functional theory calculations using spectral finite-element discretization. The proposed method treats both metallic and insulating materials in a single framework, and is applicable to both pseudopotential as well as all-electron calculations. The key ideas involved in the method include: (i) employing a higher-order spectral finite-element basis that is amenable to mesh adaption; (ii) using a Chebyshev filter to construct a subspace which is an approximation to the occupied eigenspace in a given self-consistent field iteration; (iii) using a localization procedure to construct a non-orthogonal localized basis spanning the Chebyshev filtered subspace; (iv) using a Fermi-operator expansion in terms of the subspace-projected Hamiltonian represented in the non-orthogonal localized basis to compute relevant quantities like the density matrix, electron density and band energy. We demonstrate the accuracy and efficiency of the approach on benchmark systems involving pseudopotential calculations on metallic aluminum nano-clusters up to 3430 atoms and on insulating alkane chains up to 7052 atoms, as well as all-electron calculations on silicon nano-clusters up to 3920 electrons. The benchmark studies revealed that accuracies commensurate with chemical accuracy can be obtained, and a subquadratic-scaling with system size was observed for the range of materials systems studied. In particular, for the alkane chains---close to linear-scaling is observed, whereas, for aluminum nano-clusters---the scaling is observed to be $\mathcal{O} (N^{1.46})$. For all-electron calculations on silicon nano-clusters, the scaling with the number of electrons is computed to be $\mathcal{O} (N^{1.75})$. Furthermore, significant computational savings have been realized with the proposed approach with respect to reference calculations.

We report on a quadratically convergent self-consistent field (QC-SCF) algorithm for the spin-projected unrestricted Hartree-Fock (SUHF) to mitigate the slow convergence of SUHF due to the presence of small eigenvalues in the orbital Hessian matrix. The new QC-SCF is robust and stable, allowing us to obtain the SUHF solutions very quickly. To demonstrate the applicability of the method, we present results for test systems with abundant non-dynamic correlation in comparison with the Roothaan repeated diagonalization, Pople extrapolation, and direct inversion of iterative subspace.

Hybrid system of block copolymers and homopolymers can self-assemble to form a variety of microstructures under different conditions. Based on the wormlike chain model, the self-consistent field theory for the mixing system, which is consisted of rigid-flexible diblock copolymer and the corresponding flexible homopolymer (rod-coil/coil), is presented under grand canonical ensembles. The corresponding self-consistent field algorithm is used for numerical simulation of critical micelle structure in rod-coil/coil blends. Similar structures and properties of self-assemled critical micelles in the coil-coil blends are obtained. When the volume fraction of the block increases, not only different micelle structures are formed, but also the stable shape of the micelle changes from plane to cylinder then to sphericity, and the change of the order parameters in the radial direction also demonstrates the transformation process of micelle structure.

We present an efficient quasi-Newton orbital solver optimized to reduce the number of gradient evaluations and other computational steps of comparable cost. The solver optimizes orthogonal orbitals by sequences of unitary rotations generated by the (preconditioned) limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm equipped with trust-region step restriction. The low-rank structure of the L-BFGS inverse Hessian is exploited when solving the trust-region problem. The efficiency of the proposed "Quasi-Newton Unitary Optimization with Trust-Region" (QUOTR) solver is compared to that of the standard Roothaan-Hall approach accelerated by the Direct Inversion of Iterative Subspace (DIIS), and other exact and approximate Newton solvers for mean-field (Hartree-Fock and Kohn-Sham) problems.

An advantage of modified virtual orbitals of Hartree–Fock method is discussed in the calculation of the second‐order perturbation energy. All the modified virtual orbitals can be fitted for the intermediate virtual states in the perturbation expansion, only if the molecular orbitals are expanded in terms of infinite basis functions and the set of molecular orbitals is infinite and complete. If the molecular orbitals are expanded in terms of finite basis functions, only the modified virtual orbitals with lower energies are appropriate to describe the intermediate virtual states, but the modified virtual orbitals with higher energies become inadequate. To explain the usefulness of the modified virtual orbitals, the calculation by the modified Hartree–Fock method without CI are compared with the calculation by the traditional Hartree–Fock method with complete CI.

Elementary Unitary MO Transformations and SCF Theory

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.

We report a Kohn-Sham density functional theory calculation of a system with more than 200 000 atoms and 800 000 electrons using a real-space high-order finite-difference method to investigate the electronic structure of large spherical silicon nanoclusters. Our system of choice was a 20nm large spherical nanocluster with 202 617 silicon atoms and 13 836 hydrogen atoms used to passivate the dangling surface bonds. To speed up the convergence of the eigenspace, we utilized Chebyshev-filtered subspace iteration, and for sparse matrix-vector multiplications, we used blockwise Hilbert space-filling curves, implemented in the PARSEC code. For this calculation, we also replaced our orthonormalization + Rayleigh-Ritz step with a generalized eigenvalue problem step. We utilized all of the 8192 nodes (458 752 processors) on the Frontera machine at the Texas Advanced Computing Center. We achieved two Chebyshev-filtered subspace iterations, yielding a good approximation of the electronic density of states. Our work pushes the limits on the capabilities of the current electronic structure solvers to nearly 106 electrons and demonstrates the potential of the real-space approach to efficiently parallelize large calculations on modern high-performance computing platforms.

The term biological systems may be used in reference to a wide class of polyatomic systems. They can be defined as minimal functional units which perform specific biological functions: enzymatic reactions, transport across membranes, or photosynthesis. At present, such systems as a whole are not amenable to quantum-chemistry studies because of their large size. The smallest enzymes are built of few thousands of atoms (e.g., lysozyme consists of 129 amino-acid subunits), the smallest nucleic acids are of similar size (e.g., t-RNA molecules consist of about 80 nucleotide subunits), whereas biological membranes are even larger and include different biological macromolecules embedded in a phospholipide medium. On the other hand, a common-sense definition of the term biological systems refers to any chemical molecule or molecular complex which is involved in biological or biochemical processes. The latter definition, which will be used throughout this review, covers not only complete functional units performing biological functions but also fragments of such units. Theoretical studies have provided data on properties of such fragments and have helped understanding of the biological processes at the molecular level. Depending upon the size of such fragments, they can be studied by means of various quantum-chemical methods. Molecular systems of up to a few thousands of atoms can be studied using semi-empirical methods. For the Hartree-Fock or Kohn-Sham density functional theory (DFT) calculations, the current size limit is a few hundreds of atoms. (Throughout the text, Hartree-Fock refers to ab initio Self-Consistent Field calculations using the approximation of linear combination of atomic orbitals.) When the desired accuracy requires the calculation of electron correlation at the ab initio level, only systems containing no more than few tens of atoms can be treated. Therefore, a theoretician aiming at the elucidation of biological processes by quantum-mechanical calculations faces two crucial issues. The first one is the selection of a fragment for modeling at the quantum-mechanical level. The second one is the assessment of the effects associated with parts of the system which cannot be modeled at the quantum-mechanical level. In this review, the DFT studies of biological systems are divided into two groups corresponding to different ways of addressing the second aforementioned issue.