Communication: A difference density picture for the self-consistent field ansatz Free

Self-consistent field (SCF) theory is the starting point for almost all quantum chemical techniques aimed at solving the electronic Schrödinger equation. Even at the approximate Hartree-Fock (HF) level, SCF already encapsulates much of the important chemical physics: electron pairing to form covalent bonds, and electrostatic, steric, and induction interactions between these bonds. In the context of Kohn-Sham density functional theory (KS-DFT), SCF provides a formally exact approach within a one-body picture.

Since Roothaan’s formulation of SCF,¹ many important advances have improved its computational efficiency. These include the development of efficient integral evaluation techniques in terms of Gaussian basis functions,^2–4 the introduction of direct inversion of the iterative subspace (DIIS) and related convergence stabilization recipes,^5–7 the widespread adoption of the superposition of atomic densities (SADs) guess for the starting electronic density matrix,^8,9 the maximal exploitation of spatial and/or rank sparsity in the integral tensors,^10–22 and the development of approaches to avoid or attenuate the impact of diagonalization of the Fock matrix.^23–30 

With these developments, modern SCF codes are routinely capable of treating systems with a few thousand atoms and a few tens of thousands of basis functions^31–34 (even larger computations³⁵ are possible in extreme cases). However, all SCF codes we are aware of solve for the total energy in terms of the total electronic density matrix. Because the chemically relevant quantity is normally a difference of absolute energies, this implies that enormous relative precision is required. For example, the representative case of Ribonuclease A (1860 atoms) has a total energy of approximately −5 × 10⁴ E_h which is composed of cancelling electrostatic contributions on the order of approximately ±10⁶ E_h. Since chemical applications typically require energy differences accurate to at least ∼10⁻³ (mE_h), a relative precision of better than ϵ ∼ 10⁻⁹ is required.

In this Communication, we parametrize the SCF procedure in terms of the difference between the true density (to be determined) and a SAD guess. This splits the SCF problem into two subproblems: determining the energy and external potential of the SAD guess, and subsequently converging the SCF in terms of the difference density matrix. The former problem is highly sparse and may thus be easily evaluated with high precision, whereas the latter requires much lower relative precision to obtain sub-mE_h accuracy in the final SCF energy. In the example of Ribonuclease A given above, the energy contribution due to the difference density matrix is only ∼4 × 10² E_h. We show that the lower precision requirements of our difference SCF (dSCF) approach allow us to exploit single-precision arithmetic (which can be up to 32 × more efficient than double-precision on modern stream processors, e.g., in the NVIDIA GTX 970 GPGPUs used below) and also to aggressively neglect small electron repulsion integrals. This leads to immediate speedups of up to 70% in our existing TeraChem SCF code.^36–39 

Theory: Hartree-Fock theory is usually written as a constrained functional of the density matrix D_pq, expressed in some finite spin-orbital basis {ϕ_p}

This energy represents the work involved in the coalescence of all electrons and nuclei from infinite separation. Respectively, the five terms in Eq. (1) represent (1) the energy of the external system of point charges, (2) the kinetic energy, (3) the interaction with the external system of point charges, (4) the Coulomb energy, and (5) the exchange energy. The density matrix is constrained to satisfy the usual conditions: Tr(D) = N_electron, D = D^†, D² = D.

Difference Hartree-Fock theory: As evidenced by >80 yr of application, the picture above is the most immediate and focuses attention on the total electronic density matrix D_pq as the key parameter. However, we have found it to be beneficial to cast the equations instead in terms of a nearly neutral and numerically smaller difference density matrix $DpqΔ$⁠, written relative to a superposition of atomic densities (SADs) matrix $Dpq0$⁠,

In an atomic orbital basis set, $Dpq0$ is atom-block diagonal, and thus easy to work with. We note that many authors use this SAD matrix,^8,9 or even more elaborate schemes,⁴⁰ as an SCF initial guess. However, our contribution is to recognize that writing the subsequent Hartree-Fock equations in terms of the difference density matrix (rather than the total density matrix) allows for more aggressive approximation during the SCF procedure for the difference Fock matrix. Our finding was inspired by a note in Becke and Dickson’s paper on finite-difference computation of the Coulomb operator,⁴¹ in which the authors surmised that better accuracy could possibly be obtained in their Numol program by solving for the difference Coulomb matrix rather than the full Coulomb matrix within their approximation.

The Hartree-Fock equations can now be written as a constrained functional of the difference density matrix,

The five terms carry the same meaning as in Eq. (1), but now with reference to the zeroth-order SAD system.

The first term is simply the Hartree-Fock energy of the zeroth-order external system,

The external potential is

This represents a much better starting point for electronic structure theory than the standard coalescence from infinity. Here E₀ represents the majority of the coalescence process, in which the electrons bind to nuclei to form atoms that are then positioned in the chosen molecular geometry (prior to relaxation of the electronic density). $Vpq0$ is roughly the screened nuclear potential within the molecular geometry.

The difference energy is

This represents the deformation energy inherent in the rearrangement of the electrons from atomic configurations to bonding orbitals, lone-pairs, etc.

The Fock matrix for this functional may be written as

This is formally identical to the usual Fock matrix except that it is written in terms of the external one-particle Hamiltonian $Hpq0$ and the difference Coulomb and exchange matrices, $JpqΔ$ and $KpqΔ$⁠, respectively.

The rough algorithm for dSCF is presented in Algorithm 1. Some comments on this algorithm are as follows.

• If the SAD density matrix is at least block diagonal, the external potential matrix will be a low-prefactor $O(N2)$ operation (a few times the cost of the nuclear potential) if the zero blocks in $Dpq0$ are ignored.

• If the SAD density matrix is additionally restricted to be neutral and spherical, the external potential matrix and external energy can be computed in linear scaling effort by grouping the nuclear and electronic terms in the electrostatic terms.⁴² For most neutral atoms, this potential is <10⁻¹⁰ a.u. beyond a radius of 8–10 a.u.

• The difference Coulomb and exchange matrices in Step 6 are amenable to aggressive approximations as the large, nearly cancelling positive and negative electrostatic contributions are already pre-summed in the external energy.

• The approach is additionally amenable to incremental Fock builds,¹⁰ in the usual manner. When restarting the incremental scheme (e.g., every tenth iteration), one builds $JpqΔ$ and $KpqΔ$ rather than J_pq and K_pq.

• The extension of the above approach to generalized Kohn-Sham DFT is entirely straightforward. Since the Kohn-Sham potential is nonlinear in the density, it should always be computed from the full density matrix.

Numerical testing: The dSCF approach was implemented as described above in the TeraChem GPGPU electronic structure theory package. The SAD density matrix is obtained from spherically averaged and spin-averaged unrestricted Hartree-Fock computations for each unique atom. The current algorithm for the $Vpq0$ matrix only exploits the block-diagonal structure of the $Dpq0$ matrix and thus scales as $O(N2)$⁠. The $Vpq0$ matrix is computed in mixed precision³⁶ with sufficient precision to provide for <1 μE_h error in E⁰. The dSCF iterative procedure was implemented by slightly modifying the existing SCF procedure of TeraChem and retains the ability to use dynamic precision³² and incremental Fock matrix construction. Separate user-specified precision and cutoff parameters were added to allow for different levels of approximation in $Vpq0$ and $JpqΔ/KpqΔ$⁠.

We test the dSCF approach using a molecular test set described previously.³² This test set includes closed-shell systems with up to 1860 atoms, and molecular charges ranging from −3 to +8. Note that the net charge of the system is carried into $DpqΔ$⁠, as $Dpq0$ is always taken from a neutral superposition of atomic densities. Thus, $Tr(DpqΔ)=Q$⁠, where Q is the net molecular charge. We note that using $Dpq0$ with partially charged atoms (e.g., as chosen by electronegativity equalization⁴³) often leads to a smaller E^Δ, but typically saves at most 1 or 2 SCF iterations (for tests on small molecules). For simplicity, we currently recommend using a partially charged SAD guess only in cases where the neutral guess leads to convergence problems (i.e., where the neutral SAD guess is outside the region of convergence for the DIIS procedure). Similarly, in open-shell systems, one can use a spin-polarized SAD guess to control the localization of the initial spin-density, if needed.

Total energies were converged to <1 μE_h—this required the default convergence of 10⁻⁵ maximum element in the DIIS gradient for most cases, except Ribonuclease A, which required 5 × 10⁻⁶ to achieve μE_h convergence. The 6-31G basis set is used throughout the discussion below. Geometries, raw output files, and additional numerical data are provided in the supplementary material.⁴⁴ 

Single precision SCF iterations: One of the major barriers to the development of a GPGPU SCF code is the large dichotomy between single and double precision performance on many GPGPU accelerators. Early attempts to develop an SCF iterative procedure using only single-precision J and K builds led to errors of >1 mE_h for molecules with more than a few hundred atoms, even when double precision was used for all accumulation steps. This motivated efforts by ourselves and others to develop mixed- and dynamic-precision approaches for SCF.^32,36,45 However, we find that the lower relative precision requirements in dSCF allow us to use pure single precision for J/K construction in the SCF iterations with errors of <1 mE_h for molecules with up to a few thousand atoms, as depicted in the first row of Table I.

Ultimately, our existing dynamic precision scheme provides errors of ∼1 μE_h for these systems, with essentially the same computational timings (as the vast majority of dynamic precision integrals are already computed in single precision). Thus, we continue to recommend the dynamic precision algorithm, but note that single-precision dSCF remains a compelling alternative for more exotic architectures (e.g., field-programmable gate arrays or embedded-system GPGPUs).

Sieved SCF iterations: A common tactic in integral-direct SCF codes is to truncate the pair space by a priori discarding pq pairs according to the Schwarz bound $(pq|pq)<ϵ$⁠, for a user-specified cutoff ϵ. However, inspection of the electrostatic contribution $EZ1≡12∑pqDpqVpqA+∑rs[pq|rs]Drs$ shows that very tight thresholds ϵ must be used in SCF. Otherwise, rs pairs carrying significant electronic density will be neglected without neglecting the corresponding nuclear terms in $VpqA$⁠. In this case, the molecule effectively acquires a net excess positive charge,⁴⁶ and the electrostatic energy becomes grossly inaccurate. However, with dSCF, the bulk of this electrostatic contribution is rolled into $EHF0$⁠, and the remainder is $EZ1Δ=12∑pqDpqΔVpq0+∑rs[pq|rs]DrsΔ$⁠. For dSCF, aggressive Schwarz screening is, on average, just as likely to remove rs pairs containing positive contributions as negative contributions, so the net charge of the system is less sensitive to the choice of ϵ. This is evident in our tests in TeraChem: the maximum errors in total energies are reduced by roughly 100-fold, as seen in the bottom section of Table I.

Dynamically sieved dSCF: To balance efficiency and accuracy considerations, we have found that the following “dynamically sieved dSCF” recipe works well in practice: dSCF is “burned in” to a DIIS error of 2 × 10⁻³ with single precision and a Schwarz cutoff of ϵ₁ = 10⁻⁴, at which point the Schwarz cutoff is tightened to ϵ₂ = 10⁻⁷ and dynamic precision is invoked. DIIS is restarted at this point to prevent contamination of the iterative history with numerical errors from the burn-in iterations. For the default convergence tolerance of 10⁻⁵, this yields a procedure with one or two additional SCF iterations, of which roughly half are faster burn-in iterations.

The accuracy/efficiency of dynamically sieved dSCF is demonstrated in Figure 1. The maximum error is 6.8 μE_h for Ribonuclease A, a system with 1860 atoms, 10 425 basis functions, and a total HF energy of −50 813.147 130 E_h. Note that this is a relative precision of 1.3 × 10⁻¹² in the total energy, very close to the double-precision machine epsilon. The total speedups range from 1.2 × for a medium-sized carbon nanotube system to 1.7 × for ubiquitin. Speedups in the J/K steps are generally much closer to 2 ×. For the largest systems, the invocation of dSCF exposes linear algebra steps (Fock matrix diagonalization and the DIIS orbital gradient) as rate-limiting steps. We note in passing that the J/K steps parallelize to multiple GPUs more efficiently than the linear algebra steps (even when using the MAGMA library⁴⁷ for matrix diagonalization), meaning that the linear algebra steps are rate-limiting when using dSCF with 4-8 GPUs.

Summary and conclusions: The dSCF approach described above cleanly splits the standard SCF approach into two sub-problems which are easier to solve separately. The computation of the zeroth-order energy and potential requires high precision but is atom-block sparse and can easily be made linear scaling. The subsequent convergence of the smaller difference energy relies on Fock builds of “mostly traceless” difference densities (the trace is, at most, the net charge of the system). As such, this stage of the procedure tolerates much more aggressive approximations than before, e.g., the use of single precision and/or extremely coarse Schwarz cutoffs. We have shown that sub-mE_h accuracy is attained in single-precision dSCF for up to a few thousand atoms. Moreover errors of <10 μE_h are achieved in dynamic-precision dSCF with Schwarz cutoffs of up to 10⁻⁷, for the same systems. Together, these approaches led to a dynamically sieved dSCF approach which is up to 70% faster than TeraChem’s existing dynamic-precision SCF algorithm, with errors of <10 μE_h. It should be noted that the implementation of dSCF within any existing SCF code is entirely straightforward and will produce much larger speedups in less-optimized codes (e.g., those which do not already rely heavily on dynamic precision and/or incremental Fock builds). Moreover, dSCF opens the door to apply a host of new electron repulsion integral (ERI) approximations to SCF theory: 10⁻⁹ relative error is no longer required. In the future, we will consider the use of $Dpq0$ from nearby geometries in dynamics applications (wherein the difference energies are a very favorable ≪1 E_h), the application of new tensor hypercontraction (THC)-type approximations^48,49 to efficiently resolve $FpqΔ$ and the invocation of this approach in periodic boundary conditions.

This work was supported by the National Science Foundation (No. ACI-1450179) and the Office of Naval Research (No. N00014-14-1-0590). R.M.P. acknowledges productive discussions with Chenchen Song, Professor Edward G. Hohenstein, and Professor Henrik Koch. This work was developed, in part, on the XStream computational resource, supported by the National Science Foundation Major Research Instrumentation program (No. ACI-1428930).