SCF Iteration Research Articles

Accelerating self-consistent field iterations in Kohn-Sham density functional theory using a low-rank approximation of the dielectric matrix

We present an efficient preconditioning technique for accelerating the fixed point iteration in real-space Kohn-Sham density functional theory (DFT) calculations. The preconditioner uses a low rank approximation of the dielectric matrix (LRDM) based on G\^ateaux derivatives of the residual of fixed point iteration along appropriately chosen direction functions. We develop a computationally efficient method to evaluate these G\^ateaux derivatives in conjunction with the Chebyshev filtered subspace iteration procedure, an approach widely used in large-scale Kohn-Sham DFT calculations. Further, we propose a variant of LRDM preconditioner based on adaptive accumulation of low-rank approximations from previous SCF iterations, and also extend the LRDM preconditioner to spin-polarized Kohn-Sham DFT calculations. We demonstrate the robustness and efficiency of the LRDM preconditioner against other widely used preconditioners on a range of benchmark systems with sizes ranging from $\sim$ 100-1100 atoms ($\sim$ 500--20,000 electrons). The benchmark systems include various combinations of metal-insulating-semiconducting heterogeneous material systems, nanoparticles with localized $d$ orbitals near the Fermi energy, nanofilm with metal dopants, and magnetic systems. In all benchmark systems, the LRDM preconditioner converges robustly within 20--30 iterations. In contrast, other widely used preconditioners show slow convergence in many cases, as well as divergence of the fixed point iteration in some cases. Finally, we demonstrate the computational efficiency afforded by the LRDM method, with up to 3.4$\times$ reduction in computational cost for the total ground-state calculation compared to other preconditioners.

Spin Alternation Rule in USCF for Through-Bond Magnetic Coupling─A New Look: Why and When Does It Arise and How To See It.

This work takes a new look at the spin alternation rule in unrestricted self-consistent-field (USCF) calculations in terms of structural characteristics such as periodicity, impurity location, and Coulomb exchange. For clarity, the systems considered are biradicals produced from linear conjugated hydrocarbons. Both site-parametrized Hamiltonian models for theoretical analysis and spin unrestricted density functional theory (DFT) calculations are used. Theoretical analysis leads to the following conclusions: (1) The diradical state is an excited state of a linear chain of N conjugated carbon atoms (when N is about ≤ 10). Spin alternation is a consequence of the (truncated) periodic symmetry combined with filling each closed-shell pi orbital with two electrons and each singly occupied molecular orbital (SOMO) with one electron. Spin polarization is evident in triplet (T) and broken symmetry (BS) solutions for an odd N and only in the T solution for an even N. Spin alternation is visible in the BS for an odd N and always remains muted in the calculated T. (2) For a doped chain with two radical centers, spin alternation is generally visible in the BS for an odd N. The sign of spin population on the radical centers in the BS indicates the stable spin. For radical centers separated by an odd (even) number of pz electrons, spin alternation favors T (S) state with FM (AFM) interaction. Spin oscillation remains less transparent for an even N without exchange. (3) In an unrestricted treatment with exchange, spin alternation becomes observable. Without SCF iterations, the more stable state can be identified from a clear spin oscillation in the BS. An irregular oscillation indicates a possible singlet ground state. These observations are supported by density functional calculations using the B3LYP functional and the 6-311+g(d,p) basis set on linear decapentaene diradicals with nitronyl nitroxide moieties substituted on two sets of conjugated atoms, (3,9) and (3,10). Because of the SCF procedure, one finds spin alternation in the T (BS) solution and erratic oscillation in the BS (T) solution of the 3,9 (3,10) diradical in respective equilibrium geometries. The ground state is T (S). DFT adiabatic coupling constants, SOMO energies, spin population plots, and SOMO lobe diagrams compare well with molecular electronic characteristics from theoretical analysis using Hamiltonian parameters.

Initial Maximum Overlap Method for Large Systems by the Quantum Mechanics/Extremely Localized Molecular Orbital Embedding Technique.

Quantum chemistry offers a large variety of methods to treat excited states. Many of them are based on a multireference wave function ansatz and are therefore characterized by an intrinsic complexity and high computational costs. To overcome these drawbacks and also some limitations of simpler single-reference approaches (e.g., configuration interaction singles and time-dependent density functional theory), the single-determinant Δself-consistent field-initial maximum overlap method (ΔSCF-IMOM) has been proposed. This strategy substitutes the aufbau principle with a criterion that occupies molecular orbitals at successive SCF iterations on the basis of their maximum overlap with a proper set of guess orbitals for the target excited state. In this way, it prevents the SCF process to collapse to the ground state wave function and provides excited state single Slater determinant solutions to the SCF equations. Here, we propose to extend the applicability of the IMOM to the treatment of localized excited states of large systems. To accomplish this task, we coupled it with the QM/ELMO (quantum mechanics/extremely localized molecular orbitals) strategy, a quantum mechanical embedding method in which the most chemically relevant part of the system is treated with traditional quantum chemical approaches, while the rest is described by extremely localized molecular orbitals transferred from recently constructed libraries or proper model molecules. After presenting the theoretical foundations of the new IMOM/ELMO technique, in this paper, we will show and discuss the results of preliminary test calculations carried out on both model systems (i.e., decanoic acid, decene, decapentaene, and solvated acrolein) and a system of biological interest (flavin mononucleotide in the flavodoxin protein). We observed that, for localized excited states, the new IMOM/ELMO method provides reliable results, and it reproduces the outcomes of fully IMOM calculations within the chemical accuracy threshold (i.e., 0.043 eV) by including only a limited number of atoms in the QM region. Furthermore, the first application of our embedding technique to a larger biological system gave completely plausible results in line with those obtained through more traditional quantum mechanical methods, thus opening the possibility of using the new approach in future investigations of photobiology problems.

Deep learning-based initial guess for minimum energy path calculations

An autoencoder that automatically generates an initial guess for the minimum energy pathway (MEP) calculations has been designed. Specifically, our autoencoder takes in the trajectories of molecular dynamics simulations as its input and facilitates the generation of feasible molecular coordinates. Two molecules (acetonitrile and alanine dipeptide) were tested using the nudged elastic band calculations and the results provided improvements over linear interpolation and image dependent pair potential methods in terms of the number of SCF iterations, demonstrating the utility of using an autoencoder type of an approach for MEP calculations.

GPU acceleration of all-electron electronic structure theory using localized numeric atom-centered basis functions

We present an implementation of all-electron density-functional theory for massively parallel GPU-based platforms, using localized atom-centered basis functions and real-space integration grids. Special attention is paid to domain decomposition of the problem on non-uniform grids, which enables compute- and memory-parallel execution across thousands of nodes for real-space operations, e.g. the update of the electron density, the integration of the real-space Hamiltonian matrix, and calculation of Pulay forces. To assess the performance of our GPU implementation, we performed benchmarks on three different architectures using a 103-material test set. We find that operations which rely on dense serial linear algebra show dramatic speedups from GPU acceleration: in particular, SCF iterations including force and stress calculations exhibit speedups ranging from 4.5 to 6.6. For the architectures and problem types investigated here, this translates to an expected overall speedup between 3–4 for the entire calculation (including non-GPU accelerated parts), for problems featuring several tens to hundreds of atoms. Additional calculations for a 375-atom Bi2Se3 bilayer show that the present GPU strategy scales for large-scale distributed-parallel simulations.

An Accurate Quantum-Based Approach to Explicit Solvent Effects: Interfacing the General Effective Fragment Potential Method with Ab Initio Electronic Structure Theory.

An interface between ab initio quantum mechanics (QM) methods and the general effective fragment potential (EFP2) method, QM-EFP2, is implemented in which the intermolecular interactions between a QM molecule and EFP fragments consist of Coulomb, polarization, exchange repulsion (exrep), and dispersion components. In order to ensure accuracy in the QM-EFP2 exrep interaction energy, the EFP2-EFP2 spherical Gaussian overlap (SGO) approximation is abandoned and replaced with the exact electron repulsion integrals (ERI) that are evaluated with a direct method to reduce disk usage. A Gaussian damping function for the QM-EFP2 Coulomb component damps both the EFP nuclear and electronic charges. A new overlap damping function has been implemented for the QM-EFP2 dispersion component. The current QM-EFP2 implementation has been benchmarked with the S22 and S66 data sets and demonstrates excellent agreement with symmetry-adapted perturbation theory (SAPT) for component energies and with coupled cluster theory [CCSD(T)] for the total interaction energies. Water clusters of various sizes (up to 256 water molecules) have been tested; it is shown that the QM-EFP2 method has an accuracy that is comparable to that of EFP2-EFP2. It has been shown previously that the accuracy of EFP2-EFP2 intermolecular interactions is comparable to that of second-order perturbation theory (MP2) or better. The implementation of the distributed data interface (DDI) parallelization scheme significantly improves the efficiency of QM-EFP2 calculations. The time to form the QM-EFP2 Fock operator per SCF iteration for water clusters scales linearly with the number EFP basis functions.

More robust Chebyshev filtering for SCF iteration, with applications in real-space DFT

Chebyshev filtered subspace iteration (CheFSI) is a popular and powerful solver for the eigenvalue problems embedded in electronic structure calculations. Effective use of ChFSI requires that the user choose good values for several parameters, notably the filter degree d. Finding a good value for d is typically done with user experience and some trial and error. If the user chooses a bad value for d, then convergence can be slow. For real-space computation using high-accuracy discretizations, the degree search space may be exceptionally large, making the search for a good degree difficult. We develop a method that is more robust to the consequences of bad choices of d; we also provide a heuristic-based method for choosing d automatically if the user is not able to venture a good guess. We note that only by-products of CheFSI are needed to track filter performance. We provide examples that demonstrate near-equivalence to CheFSI when d is near-optimal, and a reduction in computational resource requirements in other cases.

Optimization of convergence criteria for fragmentation methods

Fragmentation methods reproduce the result of the exponential scaling full system calculation within an error about 1kcal/mol. They serve asa robust linear scaling regime of quantum treatment for large molecular systems. As the response of fragmentation methods to convergence criteria in SCF iterations is different from that of full system calculation and the propagated convergence error is much smaller than the fragmentation error, the speed of fragmentation ab initio calculation can be further enhanced via loosening the convergence criterion in the calculation of each fragment without sacrificing the accuracy of the overall energy. Tests on single configuration calculations, geometry optimization and molecular dynamics simulation are performed.

A density functional for core-valence correlation energy.

A density functional, εCV-DFT(ρc, ρv), describing the core-valence correlation energy has been constructed as a linear combination of εLY P (corr)(ρc), εV WN5 (corr)(ρc, ρv), εPBE (corr)(ρc, ρv), εSlater (ex)(ρc, ρv), εHCTH (ex)(ρc, ρv), εHF (ex)(ρc, ρv), and FCV-DFTNi,Zi, a function of the nuclear charges. This functional, with 6 adjustable parameters, reproduces (±0.27 kcal/mol rms error) a benchmark set of 194 chemical energy changes including 9 electron affinities, 18 ionization potentials, and 167 total atomization energies covering the first- and second-rows of the periodic table. This is almost twice the rms error (±0.16 kcal/mol) obtained with CCSD(T)/MTsmall calculations, but less than half the rms error (±0.65 kcal/mol) obtained with MP2/GTlargeXP calculations, and somewhat smaller than the rms error (±0.39 kcal/mol) obtained with CCSD/MTsmall calculations. The largest positive and negative errors from εCV-DFT(ρc, ρv) were 0.88 and -0.75 kcal/mol with the set of 194 core-valence energy changes ranging from +3.76 kcal/mol for the total atomization energy of propyne to -9.05 kcal/mol for the double ionization of Mg. Evaluation of the εCV-DFT(ρc, ρv) functional requires less time than a single SCF iteration, and the accuracy is adequate for any model chemistry based on the CCSD(T) level of theory.

Extended Hückel and Slater’s rule initial guess for real space grid-based density functional theory

The Extended Hückel method is commonly used as an initial guess for Gaussian basis set quantum chemistry calculations. In this work, we combine the Extended Hückel method and Slater’s rules to form an initial guess of the electron density and molecular orbitals for all-electron real-space finite element Kohn–Sham density functional theory (FE-KSDFT). With this approach, the first SCF iteration is close in energy to the converged solution, and a factor of two speedup is obtained vs. a naïve initial electron density composed of a superposition of atomic densities. In addition, because the guessed molecular orbitals are already partially self-consistent with the guess density, we can now employ less robust, but fast, block eigensolvers such as Chebyshev-filtered subspace iteration even in the first SCF iteration to achieve 2–10× further gains vs. a conventional Jacobi–Davidson eigensolver.

A framework of constraint preserving update schemes for optimization on Stiefel manifold

This paper considers optimization problems on the Stiefel manifold $X^{\mathsf{T}}X=I_p$, where $X\in \mathbb{R}^{n \times p}$ is the variable and $I_p$ is the $p$-by-$p$ identity matrix. A framework of constraint preserving update schemes is proposed by decomposing each feasible point into the range space of $X$ and the null space of $X^{\mathsf{T}}$. While this general framework can unify many existing schemes, a new update scheme with low complexity cost is also discovered. Then we study a feasible Barzilai-Borwein-like method under the new update scheme. The global convergence of the method is established with an adaptive nonmonotone line search. The numerical tests on the nearest low-rank correlation matrix problem, the Kohn-Sham total energy minimization and a specific problem from statistics demonstrate the efficiency of the new method. In particular, the new method performs remarkably well for the nearest low-rank correlation matrix problem in terms of speed and solution quality and is considerably competitive with the widely used SCF iteration for the Kohn-Sham total energy minimization.

Self-consistent field convergence for proteins: a comparison of full and localized-molecular-orbital schemes

Proteins in the gas phase present an extreme (and unrealistic) challenge for self-consistent-field iteration schemes because their ionized groups are very strong electron donors or acceptors, depending on their formal charge. This means that gas-phase proteins have a very small band gap but that their frontier orbitals are localized compared to "normal" conjugated semiconductors. The frontier orbitals are thus likely to be separated in space so that they are close to, but not quite, orthogonal during the SCF iterations. We report full SCF calculations using the massively parallel EMPIRE code and linear scaling localized-molecular-orbital (LMO) calculations using Mopac2009. The LMO procedure can lead to artificially over-polarized wavefunctions in gas-phase proteins. The full SCF iteration procedure can be very slow to converge because many cycles are needed to overcome the over-polarization by inductive charge shifts. Example molecules have been constructed to demonstrate this behavior. The two approaches give identical results if solvent effects are included.

On a self-consistent-field-like iteration for maximizing the sum of the Rayleigh quotients

In this paper, we consider efficient methods for maximizing x⊤Bxx⊤Wx+x⊤Dx over the unit sphere, where B,D are symmetric matrices, and W is symmetric and positive definite. This problem can arise in the downlink of a multi-user MIMO system and in the sparse Fisher discriminant analysis in pattern recognition. It is already known that the problem of finding a global maximizer is closely associated with a nonlinear extreme eigenvalue problem. Rather than resorting to some general optimization methods, we introduce a self-consistent-field-like (SCF-like) iteration for directly solving the resulting nonlinear eigenvalue problem. The SCF iteration is widely used for solving the nonlinear eigenvalue problems arising in electronic structure calculations. One attractive feature of the SCF for our problem is that once it converges, the limit point not only satisfies the necessary local optimality conditions automatically, but also, and most importantly, satisfies a global optimality condition, which generally is not achievable in some optimization-based methods. The global convergence and local quadratic convergence rate are proved for certain situations. For the general case, we then discuss a trust-region SCF iteration for stabilizing the SCF iteration, which is of good global convergence behavior. Our preliminary numerical experiments show that these algorithms are more efficient than some optimization-based methods.

Sparse Projected-Gradient Method As a Linear-Scaling Low-Memory Alternative to Diagonalization in Self-Consistent Field Electronic Structure Calculations.

Large-scale electronic structure calculations usually involve huge nonlinear eigenvalue problems. A method for solving these problems without employing expensive eigenvalue decompositions of the Fock matrix is presented in this work. The sparsity of the input and output matrices is preserved at every iteration, and the memory required by the algorithm scales linearly with the number of atoms of the system. The algorithm is based on a projected gradient iteration applied to the constraint fulfillment problem. The computer time required by the algorithm also scales approximately linearly with the number of atoms (or non-null elements of the matrices), and the algorithm is faster than standard implementations of modern eigenvalue decomposition methods for sparse matrices containing more than 50 000 non-null elements. The new method reproduces the sequence of semiempirical SCF iterations obtained by standard eigenvalue decomposition algorithms to good precision.

A new theory for symmetry orbital and tensor (II)

The symmetry orbital-symmetry orbital tensor methtd is applied to the evaluation of molecular integrals (one-electron) and two-electron integrals) and the symmetry-orbital-tensor and self-consistent-field (SOT-SCF) calculations. A calculation sheme is proposed to simplify the evaluation of integrals and a key equation is derived to reduce the computation efforts in SCF iterations. According to the key equation, compared with the traditional SCF method, the computation efficiencies including CPU timing and external disk (or internal memory) requirement increase in the magnitude of the square of the order of a point group. The new SOT method is expected to be useful in the theoretical calculations of large molecular systems of high point group symmetries.

A new theory for symmetry orbital and tensor (II)——Symmetric reduction of molecular integrals and self-consistent field calculations

The symmetry orbital-symmetry orbital tensor method is applied to the evaluation of molecular integrals (one-electron and two-electron integrals) and the symmetry-orbital-tensor and self-consistent-field (SOT-SCF) calculations. A calculation scheme is proposed to simplify the evaluation of integrals and a key equation is derived to reduce the computation efforts in SCF iterations. According to the key equation, compared with the traditional SCF method, the computation efficiencies including CPU timing and external disk (or internal memory) requirement increase in the magnitude of the square of the order of a point group. The new SOT method is expected to be useful in the theoretical calculations of large molecular systems of high point group symmetries.

Self-consistent relativistic effective core potentials for transition metal atoms: Cu, Ag, and Au

A development of the RECP method for the case of transition and rare earth elements is suggested. New terms with projectors on the occupation numbers ofd andf outermost shells respectively (which can be determined in SCF iterations) are added to the standard RECP operator and the corresponding “self-consistent” RECP terms are generated for atoms Cu, Ag and Au. Significant improvement is achieved in reproduction of atomic excitation energies as compared with the conventional shape-consistent RECP calculation.

Graph theory and the PPP method

A graph theoretical formulation of the PPP method is presented. A weighted adjacency matrix of the PPP graph is given, wherein the off-diagonal elements are the bond orders. The automorphism group of the PPP graph is defined and shown to be isomorphic with the permutational subgroup of the permutation-inversion group of the molecule. It is demonstrated that the characteristic polynomial of the adjacency matrix of the PPP bond graph is invariant in every SCF iteration. It is shown that the PPP spectra discriminate isospectral graphs.

A program system for ab initio MO calculations on vector and parallel processing machines: II. SCF closed-shell and open-shell iterations

This series of three papers presents a program system for ab initio molecular orbital calculations on vector and parallel computers. Part II is devoted to SCF iterations on closed-shell and open-shell configurations starting from a file of two-electron integrals on the basis of contracted Gaussians (CGTOs). In a preliminary step, the two-electron integrals ( pq‖ rs) are reordered according to increasing values of index pq = p( p−1)/2+ q. Then, in the first SCF iteration step, a file of semi-ordered P supermatrix elements (or P and Q supermatrix elements in the open-shell case) is generated from the file of semi-ordered integrals. This file is processed at each iteration step in an efficient vector loop to generate the electron repulsion matrix. Convergence is automatically controlled through level-shifting techniques. The most time-consuming parts are the integral sorting and the generation of the P supermatrix, which are carried out only once. Subsequent SCF iteration steps respectively require 0.9 s and 1.7 s for the process of 10 7 supermatrix elements by the closed-shell and the open-shell programs on a CRAY-2 processor.

Computational chemistry with transputers: a direct SCF program

By using transputers it is possible to build up networks of parallel processors with varying topology. Due to the architecture of the processors it is appropriate to use the MIMD (multiple instruction multiple data) concept of parallel computing. The most suitable programming language is OCCAM. We investigate the use of transputer networks in computational chemistry, starting with the direct SCF method. The most time consuming step, the calculation of the two electron integrals is executed parallelly. Each node in the network calculates whole batches of integrals. The main program is written in OCCAM. For some large-scale arithmetic processes running on a single node, however, we used FORTRAN subroutines out of standard ab-initio programs to reduce the programming effort. Test calculations show, that the integral calculation step can be parallelized very efficiently. We observe a speed-up of almost 8 using eight network processors. Even in consideration of the scalar part of the SCF iteration, the speed-up is not less than 7.1.