Self-consistent Field Iteration Research Articles

Accelerating the search for global minima on potential energy surfaces using machine learning.

Controlling molecule-surface interactions is key for chemical applications ranging from catalysis to gas sensing. We present a framework for accelerating the search for the global minimum on potential surfaces, corresponding to stable adsorbate-surface structures. We present a technique using Bayesian inference that enables us to predict converged density functional theory potential energies with fewer self-consistent field iterations. We then discuss how this technique fits in with the Bayesian Active Site Calculator, which applies Bayesian optimization to the problem. We demonstrate the performance of our framework using a hematite (Fe2O3) surface and present the adsorption sites found by our global optimization method for various simple hydrocarbons on the rutile TiO2 (110) surface.

Chebyshev polynomial filtered subspace iteration in the discontinuous Galerkin method for large-scale electronic structure calculations.

The Discontinuous Galerkin (DG) electronic structure method employs an adaptive local basis (ALB) set to solve the Kohn-Sham equations of density functional theory in a discontinuous Galerkin framework. The adaptive local basis is generated on-the-fly to capture the local material physics and can systematically attain chemical accuracy with only a few tens of degrees of freedom per atom. A central issue for large-scale calculations, however, is the computation of the electron density (and subsequently, ground state properties) from the discretized Hamiltonian in an efficient and scalable manner. We show in this work how Chebyshev polynomial filtered subspace iteration (CheFSI) can be used to address this issue and push the envelope in large-scale materials simulations in a discontinuous Galerkin framework. We describe how the subspace filtering steps can be performed in an efficient and scalable manner using a two-dimensional parallelization scheme, thanks to the orthogonality of the DG basis set and block-sparse structure of the DG Hamiltonian matrix. The on-the-fly nature of the ALB functions requires additional care in carrying out the subspace iterations. We demonstrate the parallel scalability of the DG-CheFSI approach in calculations of large-scale two-dimensional graphene sheets and bulk three-dimensional lithium-ion electrolyte systems. Employing 55 296 computational cores, the time per self-consistent field iteration for a sample of the bulk 3D electrolyte containing 8586 atoms is 90 s, and the time for a graphene sheet containing 11 520 atoms is 75 s.

Numerical reduction of self‐consistent field models of macromolecular systems

AbstractWe report on the first application of model reduction techniques to the numerical simulation of self‐consistent field (SCF) models for polymer melts at equilibrium. As a complement to the classical SCF iteration, the parametric reaction‐diffusion equations for the chain propagators are treated by a reduced basis (RB) approximation from spaces of bandlimited functions. Numerical examples in a spatially one‐dimensional, periodic setting show that the approach is feasible, although not all requirements of an orthodox RB scheme are satisfied. (© 2016 Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim)

On the existence of the optimal order for wavefunction extrapolation in Born-Oppenheimer molecular dynamics.

Wavefunction extrapolation greatly reduces the number of self-consistent field (SCF) iterations and thus the overall computational cost of Born-Oppenheimer molecular dynamics (BOMD) that is based on the Kohn-Sham density functional theory. Going against the intuition that the higher order of extrapolation possesses a better accuracy, we demonstrate, from both theoretical and numerical perspectives, that the extrapolation accuracy firstly increases and then decreases with respect to the order, and an optimal extrapolation order in terms of minimal number of SCF iterations always exists. We also prove that the optimal order tends to be larger when using larger MD time steps or more strict SCF convergence criteria. By example BOMD simulations of a solid copper system, we show that the optimal extrapolation order covers a broad range when varying the MD time step or the SCF convergence criterion. Therefore, we suggest the necessity for BOMD simulation packages to open the user interface and to provide more choices on the extrapolation order. Another factor that may influence the extrapolation accuracy is the alignment scheme that eliminates the discontinuity in the wavefunctions with respect to the atomic or cell variables. We prove the equivalence between the two existing schemes, thus the implementation of either of them does not lead to essential difference in the extrapolation accuracy.

Adaptively Compressed Exchange Operator.

The Fock exchange operator plays a central role in modern quantum chemistry. The large computational cost associated with the Fock exchange operator hinders Hartree-Fock calculations and Kohn-Sham density functional theory calculations with hybrid exchange-correlation functionals, even for systems consisting of hundreds of atoms. We develop the adaptively compressed exchange operator (ACE) formulation, which greatly reduces the computational cost associated with the Fock exchange operator without loss of accuracy. The ACE formulation is not dependent on the size of the band gap, and thus can be applied to insulating and semiconducting systems, as well as metallic systems. In an iterative framework for solving Hartree-Fock-like systems, such as that observed in planewave-based methods, the ACE formulation only requires moderate modification of the code. The ACE formulation can also be advantageous for other types of basis sets, especially when the storage cost of the exchange operator is expensive. Numerical results indicate that the ACE formulation can become advantageous, even for small systems with tens of atoms. In particular, the cost of each self-consistent field iteration for the electron density in the ACE formulation is only marginally larger than that of the generalized gradient approximation (GGA) calculation, and thus offers orders-of-magnitude acceleration for Hartree-Fock-like calculations.

Accelerating Wave Function Convergence in Interactive Quantum Chemical Reactivity Studies.

The inherently high computational cost of iterative self-consistent field (SCF) methods proves to be a critical issue delaying visual and haptic feedback in real-time quantum chemistry. In this work, we introduce two schemes for SCF acceleration. They provide a guess for the initial density matrix of the SCF procedure generated by extrapolation techniques. SCF optimizations then converge in fewer iterations, which decreases the execution time of the SCF optimization procedure. To benchmark the proposed propagation schemes, we developed a test bed for performing quantum chemical calculations on sequences of molecular structures mimicking real-time quantum chemical explorations. Explorations of a set of six model reactions employing the semi-empirical methods PM6 and DFTB3 in this testing environment showed that the proposed propagation schemes achieved speedups of up to 30% as a consequence of a reduced number of SCF iterations.

Periodic Pulay method for robust and efficient convergence acceleration of self-consistent field iterations

Pulay's Direct Inversion in the Iterative Subspace (DIIS) method is one of the most widely used mixing schemes for accelerating the self-consistent solution of electronic structure problems. In this work, we propose a simple generalization of DIIS in which Pulay extrapolation is performed at periodic intervals rather than on every self-consistent field iteration, and linear mixing is performed on all other iterations. We demonstrate through numerical tests on a wide variety of materials systems in the framework of density functional theory that the proposed generalization of Pulay's method significantly improves its robustness and efficiency.

Self-consistent field for fragmented quantum mechanical model of large molecular systems.

Fragment-based linear scaling quantum chemistry methods are a promising tool for the accurate simulation of chemical and biomolecular systems. Because of the coupled inter-fragment electrostatic interactions, a dual-layer iterative scheme is often employed to compute the fragment electronic structure and the total energy. In the dual-layer scheme, the self-consistent field (SCF) of the electronic structure of a fragment must be solved first, then followed by the updating of the inter-fragment electrostatic interactions. The two steps are sequentially carried out and repeated; as such a significant total number of fragment SCF iterations is required to converge the total energy and becomes the computational bottleneck in many fragment quantum chemistry methods. To reduce the number of fragment SCF iterations and speed up the convergence of the total energy, we develop here a new SCF scheme in which the inter-fragment interactions can be updated concurrently without converging the fragment electronic structure. By constructing the global, block-wise Fock matrix and density matrix, we prove that the commutation between the two global matrices guarantees the commutation of the corresponding matrices in each fragment. Therefore, many highly efficient numerical techniques such as the direct inversion of the iterative subspace method can be employed to converge simultaneously the electronic structure of all fragments, reducing significantly the computational cost. Numerical examples for water clusters of different sizes suggest that the method shall be very useful in improving the scalability of fragment quantum chemistry methods.

DGDFT: A massively parallel method for large scale density functional theory calculations.

We describe a massively parallel implementation of the recently developed discontinuous Galerkin density functional theory (DGDFT) method, for efficient large-scale Kohn-Sham DFT based electronic structure calculations. The DGDFT method uses adaptive local basis (ALB) functions generated on-the-fly during the self-consistent field iteration to represent the solution to the Kohn-Sham equations. The use of the ALB set provides a systematic way to improve the accuracy of the approximation. By using the pole expansion and selected inversion technique to compute electron density, energy, and atomic forces, we can make the computational complexity of DGDFT scale at most quadratically with respect to the number of electrons for both insulating and metallic systems. We show that for the two-dimensional (2D) phosphorene systems studied here, using 37 basis functions per atom allows us to reach an accuracy level of 1.3 × 10(-4) Hartree/atom in terms of the error of energy and 6.2 × 10(-4) Hartree/bohr in terms of the error of atomic force, respectively. DGDFT can achieve 80% parallel efficiency on 128,000 high performance computing cores when it is used to study the electronic structure of 2D phosphorene systems with 3500-14 000 atoms. This high parallel efficiency results from a two-level parallelization scheme that we will describe in detail.

Scaling up Hartree–Fock calculations on Tianhe-2

This paper presents a new optimized and scalable code for Hartree–Fock self-consistent field iterations. Goals of the code design include scalability to large numbers of nodes, and the capability to simultaneously use CPUs and Intel Xeon Phi coprocessors. Issues we encountered as we optimized and scaled up the code on Tianhe-2 are described and addressed. A major issue is load balance, which is made challenging due to integral screening. We describe a general framework for finding a well-balanced static partitioning of the load in the presence of screening. Work stealing is used to polish the load balance. Performance results are shown on Stampede and Tianhe-2 supercomputers. Scalability is demonstrated on large simulations involving 2938 atoms and 27,394 basis functions, utilizing 8100 nodes of Tianhe-2.

On the Analysis of the Discretized Kohn--Sham Density Functional Theory

In this paper, we study a few theoretical issues in the discretized Kohn--Sham (KS) density functional theory. The equivalence between either a local or global minimizer of the KS total energy minimization problem and the solution to the KS equation is established under certain assumptions. The nonzero charge densities of a strong local minimizer are shown to be bounded from below by a positive constant uniformly. We analyze the self-consistent field (SCF) iteration by formulating the KS equation as a fixed point map with respect to the potential. The Jacobian of these fixed point maps is derived explicitly. Both global and local convergence of the simple mixing scheme can be established if the gap between the occupied states and unoccupied states is sufficiently large. This assumption can be relaxed in certain cases. Numerical experiments based on the MATLAB toolbox KSSOLV show that it holds on a few simple examples. Although our assumption on the gap is very stringent, our analysis is still valuable for a better understanding of the KS minimization problem, the KS equation, and the SCF iteration.

A Proximal Gradient Method for Ensemble Density Functional Theory

The ensemble density functional theory (E-DFT) is valuable for simulations of metallic systems due to the absence of a gap in the spectrum of the Hamiltonian matrices. Although the widely used self-consistent field (SCF) iteration method can be extended to solve the minimization of the total energy functional with respect to orthogonality constraints, there is no theoretical guarantee on the convergence of these algorithms. In this paper, we consider an equivalent model with a single variable and a single spherical constraint by eliminating the dependence on the fractional occupancies. A proximal gradient method is developed by keeping the entropy term but linearizing all other terms in the total energy functional. Convergence to the stationary point is established. Numerical results using the MATLAB toolbox KSSOLV can outperform SCF consistently on many metallic systems.

Variational, Self-Consistent Implementation of the Perdew-Zunger Self-Interaction Correction with Complex Optimal Orbitals.

A variational, self-consistent implementation of the Perdew-Zunger self-interaction correction (PZ-SIC), based on a unified Hamiltonian and complex optimal orbitals, is presented for finite systems and atom-centered basis sets. A simplifying approximation allowing the use of real canonical orbitals is proposed. The algorithm is based on two-step self-consistent field iterations, where the updates of the canonical orbitals and the optimal orbitals are done separately. Calculations of the energy of atoms ranging from H to Ar are presented, using various generalized gradient functionals (PBE, APBE, PBEsol) and a meta-generalized gradient functional (TPSS). While the energy of atoms is poorly described by PBEsol, which is a functional optimized to reproduce properties of solids, the PZ-SIC brings the calculations into good agreement with the best ab initio estimates. The importance of using complex optimal orbitals becomes particularly clear in calculations using the TPSS functional, where the original functional gives good results while the application of PZ-SIC with real orbitals gives highly inaccurate results. With complex optimal orbitals, PZ-SIC slightly improves the accuracy of the TPSS functional. The charge localization problem that plagues Kohn-Sham DFT functionals, including hybrid functionals, is illustrated by calculations on the CH3 + F(-) complex, where even PBEsol with PZ-SIC is found to give estimates of both energy and charge with accuracy comparable to that of coupled cluster calculations.

Multiple environment single system quantum mechanical/molecular mechanical (MESS-QM/MM) calculations. 1. Estimation of polarization energies.

In combined quantum mechanical/molecular mechanical (QM/MM) free energy calculations, it is often advantageous to have a frozen geometry for the quantum mechanical (QM) region. For such multiple-environment single-system (MESS) cases, two schemes are proposed here for estimating the polarization energy: the first scheme, termed MESS-E, involves a Roothaan step extrapolation of the self-consistent field (SCF) energy; whereas the other scheme, termed MESS-H, employs a Newton–Raphson correction using an approximate inverse electronic Hessian of the QM region (which is constructed only once). Both schemes are extremely efficient, because the expensive Fock updates and SCF iterations in standard QM/MM calculations are completely avoided at each configuration. They produce reasonably accurate QM/MM polarization energies: MESS-E can predict the polarization energy within 0.25 kcal/mol in terms of the mean signed error for two of our test cases, solvated methanol and solvated β-alanine, using the M06-2X or ωB97X-D functionals; MESS-H can reproduce the polarization energy within 0.2 kcal/mol for these two cases and for the oxyluciferin–luciferase complex, if the approximate inverse electronic Hessians are constructed with sufficient accuracy.

Subquadratic-scaling subspace projection method for large-scale Kohn-Sham density functional theory calculations using spectral finite-element discretization

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.

Maximization of the sum of the trace ratio on the Stiefel manifold, II: Computation

The necessary condition established in Part I of this paper for the global maximizers of the maximization problem $$\mathop {\max }\limits_V \left\{ {\frac{{tr(V^ \top AV)}} {{tr(V^ \top BV)}} + tr(V^ \top CV)} \right\} $$ over the Stiefel manifold \(\{ V \in \mathbb{R}^{m \times \ell } |V^ \top V = I_\ell \} (\ell < m)\), naturally leads to a self-consistent-field (SCF) iteration for computing a maximizer. In this part, we analyze the global and local convergence of the SCF iteration, and show that the necessary condition for the global maximizers is fulfilled at any convergent point of the sequences of approximations generated by the SCF iteration. This is one of the advantages of the SCF iteration over optimization-based methods. Preliminary numerical tests are reported and show that the SCF iteration is very efficient by comparing with some manifold-based optimization methods.

Acceleration of self-consistent field convergence in ab initio molecular dynamics simulation with multiconfigurational wave function.

The Lagrange interpolation of molecular orbital (LIMO) method, which reduces the number of self-consistent field iterations in ab initio molecular dynamics simulations with the Hartree-Fock method and the Kohn-Sham density functional theories, is extended to the theory of multiconfigurational wave functions. We examine two types of treatments for the active orbitals that are partially occupied. The first treatment, as denoted by LIMO(C), is a simple application of the conventional LIMO method to the union of the inactive core and the active orbitals. The second, as denoted by LIMO(S), separately treats the inactive core and the active orbitals. Numerical tests to compare the two treatments clarify that LIMO(S) is superior to LIMO(C). Further applications of LIMO(S) to various systems demonstrate its effectiveness and robustness.

On the Convergence of the Self-Consistent Field Iteration in Kohn--Sham Density Functional Theory

It is well known that the self-consistent field (SCF) iteration for solving the Kohn--Sham (KS) equation often fails to converge, yet there is no clear explanation. In this paper, we investigate the SCF iteration from the perspective of minimizing the corresponding KS total energy functional. By analyzing the second-order Taylor expansion of the KS total energy functional and estimating the relationship between the Hamiltonian and the part of the Hessian which is not used in the SCF iteration, we are able to prove global convergence from an arbitrary initial point and local linear convergence from an initial point sufficiently close to the solution of the KS equation under the assumptions that the gap between the occupied states and unoccupied states is sufficiently large and the second-order derivatives of the exchange correlation functional are uniformly bounded from above. Although these conditions are very stringent and are almost never satisfied in reality, our analysis is interesting in the sense th...

A third-generation density-functional-theory-based method for calculating canonical molecular orbitals of large molecules

We used grid-free modified Cholesky decomposition (CD) to develop a density-functional-theory (DFT)-based method for calculating the canonical molecular orbitals (CMOs) of large molecules. Our method can be used to calculate standard CMOs, analytically compute exchange-correlation terms, and maximise the capacity of next-generation supercomputers. Cholesky vectors were first analytically downscaled using low-rank pivoted CD and CD with adaptive metric (CDAM). The obtained Cholesky vectors were distributed and stored on each computer node in a parallel computer, and the Coulomb, Fock exchange, and pure exchange-correlation terms were calculated by multiplying the Cholesky vectors without evaluating molecular integrals in self-consistent field iterations. Our method enables DFT and massively distributed memory parallel computers to be used in order to very efficiently calculate the CMOs of large molecules.

Gradient Type Optimization Methods For Electronic Structure Calculations

The density functional theory (DFT) in electronic structure calculations can be formulated as either a nonlinear eigenvalue or a direct minimization problem. The most widely used approach for solving the former is the so-called self-consistent field (SCF) iteration. A common observation is that the convergence of SCF is not clear theoretically, while approaches with convergence guarantees for solving the latter are often not competitive with SCF numerically. In this paper, we study gradient type methods for solving the direct minimization problem by constructing new iterations along the gradient on the Stiefel manifold. Global convergence rate (i.e., convergence to a stationary point from any initial solution) as well as local convergence rate can be obtained from the standard theory for gradient methods on manifold directly. A major computational advantage is that the computation of linear eigenvalue problems is no longer needed. The main costs of our approaches arise from the assembling of the total ene...