We formulate self-consistent field (SCF) theory in terms of an interaction picture where the working variable is the difference density matrix between the true system and a corresponding superposition of atomic densities. As the difference density matrix directly represents the electronic deformations inherent in chemical bonding, this “difference self-consistent field (dSCF)” picture provides a number of significant conceptual and computational advantages. We show that this allows for a stable and efficient dSCF iterative procedure with wholly single-precision Coulomb and exchange matrix builds. We also show that the dSCF iterative procedure can be performed with aggressive screening of the pair space. These approximations are tested and found to be accurate for systems with up to 1860 atoms and >10 000 basis functions, providing for immediate overall speedups of up to 70% in the heavily optimized TeraChem SCF implementation.

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,1 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 functions31–34 (even larger computations35 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 × 104 Eh which is composed of cancelling electrostatic contributions on the order of approximately ±106 Eh. Since chemical applications typically require energy differences accurate to at least ∼10−3 (mEh), a relative precision of better than ϵ ∼ 10−9 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-mEh 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 × 102 Eh. 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 Dpq, expressed in some finite spin-orbital basis {ϕp}

EHF(Dpq)EA+pqDpqTpq+DpqVpqA+12pqrsDpq[pq|rs]DrsDpq[pr|qs]Drs.
(1)

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) = Nelectron, D = D, D2 = 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 Dpq 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,

DpqDpq0+DpqΔ.
(2)

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,40 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,41 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,

EHF(DpqΔ)E0+pqDpqΔTpq+DpqΔVpq0+12pqrsDpqΔ[pq|rs]DrsΔDpqΔ[pr|qs]DrsΔ.
(3)

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,

E0EA+pqDpq0Tpq+Dpq0VpqA+12pqrsDpq0[pq|rs]Drs0Dpq0[pr|qs]Drs0.
(4)

The external potential is

Vpq0VpqA+rs[pq|rs]Drs0[pr|qs]Drs0.
(5)

This represents a much better starting point for electronic structure theory than the standard coalescence from infinity. Here E0 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

EΔ(DpqΔ)pqDpqΔTpq+DpqΔVpq0+12pqrsDpqΔ[pq|rs]DrsΔDpqΔ[pr|qs]DrsΔ.
(6)

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

FpqEΔDpqΔ=Tpq+Vpq0+rs[pq|rs]DrsΔ[pr|qs]DrsΔHpq0+JpqΔKpqΔ.
(7)

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.

ALGORITHM 1.

Generic dSCF Algorithm.

 
1: proceduredSCF (Molecule, Basis Set)  
2: Obtain Dpq0 ⊳ SAD Guess 
3: Compute External Potential: ⊳ Sparse 
Vpq0=VpqA+[pq|rs]Drs0[pr|qs]Drs0 
4: Compute External Energy: ⊳ Large 
E0EA+Dpq0Tpq+Dpq0VpqA 
+12Dpq0[pq|rs]Drs012Dpq0[pr|qs]Drs0 
5: Compute External Core Hamiltonian: 
Hpq0=Tpq+Vpq0 
6: Set DpqΔ=0 ⊳ Optionally, use other guess 
7: while Not Converged do  
8: Compute Difference J/K Matrices: ⊳ Dense 
JpqΔ=[pq|rs]DrsΔ,KpqΔ=[pr|qs]DrsΔ 
9: Compute Fock Matrix:  
Fpq=Hpq0+JpqΔKpqΔ 
10: Compute Difference Energy: ⊳ Small 
EΔ=12DpqΔHpq0+Fpq 
11: Diagonalize Fock Matrix:  
FpqCqi=SpqCqiϵi 
12: Compute Difference Density Matrix: 
DpqΔ=CpiCqiDpq0 
13: end while  
14: end procedure  
 
1: proceduredSCF (Molecule, Basis Set)  
2: Obtain Dpq0 ⊳ SAD Guess 
3: Compute External Potential: ⊳ Sparse 
Vpq0=VpqA+[pq|rs]Drs0[pr|qs]Drs0 
4: Compute External Energy: ⊳ Large 
E0EA+Dpq0Tpq+Dpq0VpqA 
+12Dpq0[pq|rs]Drs012Dpq0[pr|qs]Drs0 
5: Compute External Core Hamiltonian: 
Hpq0=Tpq+Vpq0 
6: Set DpqΔ=0 ⊳ Optionally, use other guess 
7: while Not Converged do  
8: Compute Difference J/K Matrices: ⊳ Dense 
JpqΔ=[pq|rs]DrsΔ,KpqΔ=[pr|qs]DrsΔ 
9: Compute Fock Matrix:  
Fpq=Hpq0+JpqΔKpqΔ 
10: Compute Difference Energy: ⊳ Small 
EΔ=12DpqΔHpq0+Fpq 
11: Diagonalize Fock Matrix:  
FpqCqi=SpqCqiϵi 
12: Compute Difference Density Matrix: 
DpqΔ=CpiCqiDpq0 
13: end while  
14: end procedure  

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.42 For most neutral atoms, this potential is <10−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,10 in the usual manner. When restarting the incremental scheme (e.g., every tenth iteration), one builds JpqΔ and KpqΔ rather than Jpq and Kpq.

  • 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 precision36 with sufficient precision to provide for <1 μEh error in E0. The dSCF iterative procedure was implemented by slightly modifying the existing SCF procedure of TeraChem and retains the ability to use dynamic precision32 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.32 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 equalization43) 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 μEh—this required the default convergence of 10−5 maximum element in the DIIS gradient for most cases, except Ribonuclease A, which required 5 × 10−6 to achieve μEh 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.44 

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 mEh 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 mEh for molecules with up to a few thousand atoms, as depicted in the first row of Table I.

TABLE I.

Maximum absolute errors (Eh) in total Hartree-Fock energies for SCF and dSCF approaches with various approximations over the SCF test set. RHF/6-31G converged to <1 μEh precision.

ApproximationSCF errordSCF error
Precision
Single precision 2.6 × 10−2 5.7 × 10−4 
Dynamic precision 2.0 × 10−6 1.0 × 10−6 
Schwarz cutoffa 
ϵ = 10−4 9.1 × 10+1 3.4 × 10−2 
ϵ = 10−5 1.4 × 10+1 1.9 × 10−3 
ϵ = 10−6 1.8 × 10−2 7.0 × 10−5 
ϵ = 10−7 8.4 × 10−4 7.2 × 10−6 
ApproximationSCF errordSCF error
Precision
Single precision 2.6 × 10−2 5.7 × 10−4 
Dynamic precision 2.0 × 10−6 1.0 × 10−6 
Schwarz cutoffa 
ϵ = 10−4 9.1 × 10+1 3.4 × 10−2 
ϵ = 10−5 1.4 × 10+1 1.9 × 10−3 
ϵ = 10−6 1.8 × 10−2 7.0 × 10−5 
ϵ = 10−7 8.4 × 10−4 7.2 × 10−6 
a

The default Schwarz cutoff in TeraChem is ϵ = 10−13, which is used throughout for comparison.

Ultimately, our existing dynamic precision scheme provides errors of ∼1 μEh 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 EZ112pqDpqVpqA+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,46 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Δ=12pqDpqΔ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−3 with single precision and a Schwarz cutoff of ϵ1 = 10−4, at which point the Schwarz cutoff is tightened to ϵ2 = 10−7 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−5, 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 μEh for Ribonuclease A, a system with 1860 atoms, 10 425 basis functions, and a total HF energy of −50 813.147 130 Eh. Note that this is a relative precision of 1.3 × 10−12 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 library47 for matrix diagonalization), meaning that the linear algebra steps are rate-limiting when using dSCF with 4-8 GPUs.

FIG. 1.

Demonstration of the error/speedup profile for dynamically sieved dHF (dHF-dyn2) vs. the default dynamic precision HF (HF-dyn) algorithm in TeraChem, for the five largest members of the SCF test set. Timings were performed using a single NVIDIA GeForce GTX 970 GPGPU.

FIG. 1.

Demonstration of the error/speedup profile for dynamically sieved dHF (dHF-dyn2) vs. the default dynamic precision HF (HF-dyn) algorithm in TeraChem, for the five largest members of the SCF test set. Timings were performed using a single NVIDIA GeForce GTX 970 GPGPU.

Close modal

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-mEh accuracy is attained in single-precision dSCF for up to a few thousand atoms. Moreover errors of <10 μEh are achieved in dynamic-precision dSCF with Schwarz cutoffs of up to 10−7, 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 μEh. 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−9 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 Eh), the application of new tensor hypercontraction (THC)-type approximations48,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).

1.
C. C. J.
Roothaan
,
Rev. Mod. Phys.
23
,
69
(
1951
).
2.
L. E.
McMurchie
and
E. R.
Davidson
,
J. Comput. Phys.
26
,
218
(
1978
).
3.
J.
Rys
,
M.
Dupuis
, and
H. F.
King
,
J. Comput. Chem.
4
,
154
(
1983
).
4.
S.
Obara
and
A.
Saika
,
J. Chem. Phys.
84
,
3963
(
1986
).
5.
6.
K. N.
Kudin
,
G. E.
Scuseria
, and
E.
Cancs
,
J. Chem. Phys.
116
,
8255
(
2002
).
7.
X.
Hu
and
W.
Yang
,
J. Chem. Phys.
132
,
054109
(
2010
).
8.
J.
Almlöf
,
K.
Faegri
,
M.
Feyereisen
, and
K.
Korsell
,
DISCO, a DirectSCF and MP2 Program
(
University of Minnesota
,
Minneapolis, MN
,
1982
).
9.
J. H.
Van Lenthe
,
R.
Zwaans
,
H. J. J.
Van Dam
, and
M. F.
Guest
,
J. Comput. Chem.
27
,
926
(
2006
).
10.
J.
Almlöf
,
K.
Faegri
, and
K.
Korsell
,
J. Comput. Chem.
3
,
385
(
1982
).
11.
R. A.
Friesner
,
J. Chem. Phys.
86
,
3522
(
1987
).
12.
M.
Häser
and
R.
Ahlrichs
,
J. Comput. Chem.
10
,
104
(
1989
).
13.
O.
Vahtras
,
J.
Almlöf
, and
M. W.
Feyereisen
,
Chem. Phys. Lett.
213
,
514
(
1993
).
14.
C. A.
White
,
B. G.
Johnson
,
P. M.
Gill
, and
M.
Head-Gordon
,
Chem. Phys. Lett.
230
,
8
(
1994
).
15.
M. C.
Strain
,
G. E.
Scuseria
, and
M. J.
Frisch
,
Science
271
,
51
(
1996
).
16.
C. A.
White
and
M.
HeadGordon
,
J. Chem. Phys.
104
,
2620
(
1996
).
17.
E.
Schwegler
,
M.
Challacombe
, and
M.
Head-Gordon
,
J. Chem. Phys.
106
,
9708
(
1997
).
18.
R.
Polly
,
H.-J.
Werner
,
F. R.
Manby
, and
P. J.
Knowles
,
Mol. Phys.
102
,
2311
(
2004
).
19.
L.
Füsti-Molnár
and
P.
Pulay
,
J. Chem. Phys.
117
,
7827
(
2002
).
20.
H.
Koch
,
A. S.
de Meras
, and
T. B.
Pedersen
,
J. Chem. Phys.
118
,
9481
(
2003
).
21.
D. S.
Lambrecht
and
C.
Ochsenfeld
,
J. Chem. Phys.
123
,
184101
(
2005
).
22.
F.
Neese
,
F.
Wennmohs
,
A.
Hansen
, and
U.
Becker
,
Chem. Phys.
356
,
98
(
2009
).
23.
R.
McWeeny
,
Proc. R. Soc. London, Ser. A
235
,
496
(
1956
).
24.
X.-P.
Li
,
R. W.
Nunes
, and
D.
Vanderbilt
,
Phys. Rev. B
47
,
10891
(
1993
).
25.
26.
J. M.
Millam
and
G. E.
Scuseria
,
J. Chem. Phys.
106
,
5569
(
1997
).
27.
A. H. R.
Palser
and
D. E.
Manolopoulos
,
Phys. Rev. B
58
,
12704
(
1998
).
28.
M.
Challacombe
,
J. Chem. Phys.
110
,
2332
(
1999
).
29.
A. D.
Daniels
and
G. E.
Scuseria
,
J. Chem. Phys.
110
,
1321
(
1999
).
30.
E. H.
Rubensson
,
E.
Rudberg
, and
P.
Sałek
,
J. Chem. Phys.
128
,
074106
(
2008
).
31.
G. E.
Scuseria
,
J. Phys. Chem. A
103
,
4782
(
1999
).
32.
N.
Luehr
,
I. S.
Ufimtsev
, and
T. J.
Martínez
,
J. Chem. Theory Comput.
7
,
949
(
2011
).
33.
F.
Neese
,
WIREs: Comput. Mol. Sci.
2
,
73
(
2012
).
34.
J.
Kussmann
,
M.
Beer
, and
C.
Ochsenfeld
,
WIREs: Comput. Mol. Sci.
3
,
614
(
2013
).
35.
E.
Rudberg
,
E. H.
Rubensson
, and
P.
Saek
,
J. Chem. Theory Comput.
7
,
340
(
2011
).
36.
I. S.
Ufimtsev
and
T. J.
Martínez
,
J. Chem. Theory Comput.
4
,
222
(
2008
).
37.
I. S.
Ufimtsev
and
T. J.
Martínez
,
J. Chem. Theory Comput.
5
,
1004
(
2009
).
38.
I. S.
Ufimtsev
and
T. J.
Martínez
,
J. Chem. Theory Comput.
5
,
2619
(
2009
).
39.
A. V.
Titov
,
I. S.
Ufimtsev
,
N.
Luehr
, and
T. J.
Martínez
,
J. Chem. Theory Comput.
9
,
213
(
2013
).
40.
B.
Jansik
 et al,
Phys. Chem. Chem. Phys.
11
,
5805
(
2009
).
41.
A. D.
Becke
and
R. M.
Dickson
,
J. Chem. Phys.
89
,
2993
(
1988
).
42.
J. M.
Soler
 et al,
J. Phys.: Condens. Matter
14
,
2745
(
2002
).
43.
W. J.
Mortier
,
S. K.
Ghosh
, and
S.
Shankar
,
J. Am. Chem. Soc.
108
,
4315
(
1986
).
44.
See supplementary material at http://dx.doi.org/10.1063/1.4945277 for test set geometries, raw input/output file pairs, and additional analysis.
45.
K.
Yasuda
,
J. Chem. Theory Comput.
4
,
1230
(
2008
).
46.
D. L.
Wilhite
and
R. N.
Euwema
,
J. Chem. Phys.
61
,
375
(
1974
).
47.
S.
Tomov
,
J.
Dongarra
, and
M.
Baboulin
,
Parallel Comput.
36
,
232
(
2010
).
48.
E. G.
Hohenstein
,
R. M.
Parrish
, and
T. J.
Martínez
,
J. Chem. Phys.
137
,
044103
(
2012
).
49.
R. M.
Parrish
,
E. G.
Hohenstein
,
T. J.
Martínez
, and
C. D.
Sherrill
,
J. Chem. Phys.
137
,
224106
(
2012
).

Supplementary Material