Approaching the basis set limit for DFT calculations using an environment-adapted minimal basis with perturbation theory: Formulation, proof of concept, and a pilot implementation

Kohn-Sham density functional theory^1–3 (KS-DFT) has become the most widespread electronic-structure method because of its reasonable balance between accuracy and computational cost. Functionals using the generalized gradient approximation (GGA)^4,5 are usually regarded as the simplest that can give acceptable accuracy for chemistry. To overcome the plague of self-interaction error, new variables have been introduced, leading to meta-GGA,^6–8 global hybrid (GH),⁹ and range-separated hybrid (RSH)^10,11 functionals. At the same time, a variety of models have been developed to account for van der Waals (vdW) interactions within DFT,^12,13 including the empirical DFT-D methods^14–16 and nonlocal correlation (NLC) functionals (e.g., vdw-DF-10,¹⁷ VV10¹⁸). Most recently, two combinatorially designed functionals were developed: ωB97X-V¹⁹ (RSH+VV10) and B97M-V²⁰ (meta-GGA+VV10), which demonstrated impressive accuracy on both thermochemistry (TC) and non-covalent interactions (NC), with an accessible complete basis set (CBS) limit, and low grid sensitivity.

With finite atomic orbital (AO) basis sets,²¹ one prerequisite for attaining such accuracy is to approach the CBS limit. This issue has been carefully investigated,^22–25 but is often neglected in practical applications, as exemplified by the prevalence of the “B3LYP/6-31G(d)” model chemistry. A basis set of at least triple-ζ and preferably quadruple-ζ quality is often required by hybrid functionals (e.g., B3LYP^9,26,27) to obtain adequately converged thermochemistry results. Even for the semi-local B97M-V functional, the acceptable alternatives to aug-cc-pVQZ (aQZ)^28,29 (which almost represents the CBS limit) are still of least triple-ζ quality. Turning to the evaluation of NCs, a similar study on the A24³⁰ and S66^31,32 complexes indicates that augmented triple-ζ basis sets (e.g., aug-cc-pVTZ (aTZ), def2-TZVPD³³) are in general required by B97M-V to properly converge the binding energies. Their double-ζ counterparts (e.g., aug-cc-pVDZ (aDZ), def2-SVPD) should only be carefully used with counterpoise (CP) corrections.³⁴ For ωB97X-V, the requirement on basis set qualities might be even higher due to the slower basis set convergence of functionals that contain exact exchange.

Figure 1 demonstrates the basis set convergence of several modern density functionals in terms of their root-mean-square deviations (RMSDs) for the G2 set³⁵ (atomization energies of 148 neutral molecules, whose reference values are taken from Ref. 36). Apart from B97M-V, using aQZ instead of aTZ for the other four functionals reduces their RMSDs by 0.8–1.4 kcal/mol, including the semi-local B97-D¹⁵ functional (GGA). Using aDZ yields very poor accuracy (usually over 10 kcal/mol RMSDs) for all these functionals, which defeats the purpose of using state-of-the-art density functionals. One way to tackle this problem is by directly training a functional in a small basis, such as the EDF1 functional,³⁷ which was parameterized at the 6-31+G(d) level. By relying on cancellation between functional error and basis set error, the transferability of these methods are often limited, and further empirical corrections seem necessary to achieve adequate accuracy for relative energies.³⁸ 

Each self-consistent field (SCF) cycle of a KS-DFT calculation involves two computationally demanding steps: (1) Fock matrix construction with a given density and (2) Fock matrix diagonalization to update the density. For fixed system size, the computational cost of the Fock build scales as $O(n4)$ with respect to the basis size (n) when conventional AO algorithms are used, and the cost of the diagonalization step scales as $O(n3)$⁠. This steep cost increase inhibits large basis sets (e.g., those of quadruple-ζ size) from being routinely employed in DFT calculations. It should be noted that the scaling of cost vs. basis size is largely independent of the development of linear scaling (with system size) Fock matrix build algorithms^40–47 and many diagonalization replacements.^48,49 Moreover, near-complete basis sets are not favored by linear-scaling algorithms, especially when diffuse functions are included, since matrix element sparsity is diminished and the overlap matrix starts to be ill-conditioned (note that optimization of Gaussian basis sets up to triple-ζ quality with much reduced condition numbers while maintaining condensed phase accuracy has been reported⁵⁰), which in turn destroys the sparsity of the density matrix.^51,52

One successful strategy to make large basis KS-DFT calculations more tractable is to compute the full Coulomb (J) and exchange (K) matrices more efficiently by approximating two-electron repulsion integrals (ERIs) with the aid of auxiliary basis functions or grid points. The resolution-of-the-identity (RI) method^53–55 expands the product of AO function pairs with a preoptimized auxiliary basis. RI algorithms do not improve the system-size scaling unless local fit regions are applied,^56,57 but they reduce the basis set size scaling from $O(n4)$ to $O(n3)$⁠. Therefore, state-of-the-art RI algorithms (e.g., MARI-J,⁵⁸ PARI-K,⁵⁹ and occ-RI-K⁶⁰) can speed up Fock matrix constructions in large basis sets significantly for small- to medium-sized systems, while retaining numerical accuracy. The diagonalization step is unaffected.

A second successful approach to accelerating large basis calculations is to perform the iterative SCF procedures in a primary (small) basis and then approximate the secondary (target) basis results by utilizing perturbation theory. This idea was introduced for post-SCF methods (e.g., MP2),^61,62 and was then developed for SCF methods, including the dual-basis SCF (DB-SCF) method developed by Head-Gordon and co-workers,^63–67 and the “Hartree-Fock/Density Functional Perturbation Corrections” (HFPC/DFPC) scheme proposed by Deng and Gill.^68–70 With a careful choice of primary/secondary basis set pairing, these methods can provide satisfactory accuracy for both TC^63,69,70 and NC⁶⁷ with significantly reduced computational costs (roughly 10 times faster), although system-size scaling remains unchanged. One limitation is the need to develop and validate the basis set pairings,⁶⁷ which determines accuracy and speedup. Furthermore, as the secondary basis approaches the CBS limit, the size of the primary basis needed to achieve a given accuracy also increases: for instance the optimized primary basis for cc-pVQZ is roughly of cc-pVTZ size.⁶⁴ 

A related approach is the use of small adaptive basis sets. The idea of encoding chemical environment information into atomic/quasiatomic basis functions to understand chemical bonding dates back to early tools,^71–76 as well as some more contemporary methods.^77–82 Apart from interpretive purposes, the merits of utilizing small adaptive bases in KS-DFT calculations have been recognized with the development of fast (especially linear-scaling) SCF algorithms, leading to renewed interest in the concept of “polarized atomic orbitals” (PAOs), first put forward by Adams in the 1960s.^83,84 These adaptive sets usually have very tiny (often minimal) spans, which leads to vastly fewer variational degrees of freedom. In addition, an adaptive basis constructed with spatial confinement contributes to a well-conditioned overlap matrix, which is a property favored by $O(N)$ scaling methods.

The Adams PAO scheme treats atoms in a molecule as fragments and solves projected equations self-consistently on each of them, which is similar, in spirit, to projected SCF methods using fragment-localized, non-orthogonal molecular orbitals (MOs) to evaluate intermolecular interactions (SCF-MI).^85–87 In practice this scheme only works for weakly interacting atoms (e.g., rare-gas clusters) or ionic compounds (e.g., LiH, NaCl).⁸⁸ Later, the PAO approach was recast to form a minimal atom-centered adaptive basis as an atom-blocked contraction of the secondary basis functions on each atom.⁸⁹ The molecular energy is minimized simultaneously with respect to the atom-blocked contraction coefficients and the density matrix in the adaptive basis.⁸⁹ The PAO-SCF energy can be improved using perturbation theory,⁹⁰ similar to the dual basis approaches discussed above. The minimal rank of the PAO basis and its atomic locality makes it promising for linear-scaling algorithms,⁵² but the “double” optimization problem is challenging and often causes convergence problems.

Significant progress on tractable adaptive basis schemes for KS-DFT has been made in the condensed matter physics community. Similar to the aforementioned PAO approach, Ozaki and Kino^91,92 and others^93,94 used numerical solutions to the atomic Kohn-Sham equations as the secondary basis, and a scheme resembling geometry optimization to obtain the adaptive basis. The CONQUEST program⁹⁵ forms local “support functions” (an adaptive basis) from either functions akin to plane waves⁹⁶ or pseudo-atomic orbitals.⁹⁷ The ONETEP package^98,99 forms non-orthogonal generalized Wannier functions (NGWFs)¹⁰⁰ as the environment-adapted basis, which is a linear combination of periodic sinc functions confined in an atom-centered sphere of fixed radius. The NGWFs are efficiently optimized via a preconditioned conjugate-gradient algorithm.¹⁰¹ 

Recently, adaptive basis schemes that do not require the global Hamiltonian or density matrix have been presented. The localized filter diagonalization (LFD) method builds an adaptive basis on-the-fly by contracting the atomic Gaussian functions within a local region, with contraction coefficients determined by diagonalizing a block of the Hamiltonian matrix corresponding to that region.^102,103 This algorithm has also been used to construct multisite local support functions,¹⁰⁴ and the general philosophy has been extended by Lin et al.,¹⁰⁵ including another model with more rigorous optimization.¹⁰⁶ While clearly promising, to our knowledge, the accuracy and performance of these methods on chemical systems have not been systematically assessed yet.

In the present work, we propose an inexpensive version of the PAO method (Sec. II). Instead of energy-optimizing the adaptive basis and density simultaneously, an inexpensive converged SCF solution (density matrix) computed in a projected reference basis (PRB) is utilized as a reference (Sec. II A). Based on this reference, an atom-centered minimal adaptive basis (MAB) is found by minimizing a judiciously chosen surrogate function (Sec. II B), which only involves computationally inexpensive steps. The converged MAB is then used as the basis set for another SCF calculation, which requires small computational effort as well while providing comparable accuracy to PAO-SCF. Perturbation corrections (PCs) can be applied to the MAB-SCF energy for obtaining the desired accuracy (Sec. II E). The overall MAB-SCF (PC) procedure is illustrated in Figure 2. Details about the pilot implementation of this scheme and proof-of-concept calculations are summarized in Sec. III. As an approximate SCF method, its accuracy is assessed on a broad range of TC and NC datasets that is presented in Secs. IV and V.

The notation used throughout this paper is as follows: |ω〉: generic atomic basis functions; |ψ〉: generic molecular orbitals; capital Roman indices X, Y, …: atomic centers; lowercase Greek letters μ, ν, λ, …: secondary (large) AO basis indices; α, β, γ, …: primary (PRB or MAB) AO basis indices; lowercase Romans i, j, k, …: occupied MO indices; a, b, c, …: virtual MO indices; p, q, r, …: generic MO indices. For introducing the MAB optimization scheme, i, j, … are also employed to denote the vectors retained in the MAB subspace, a, b, … for the vectors in MAB’s complementary subspace, and p, q, … for the generic ones, which is analogous to the partitioning of MO space in SCF. The different basis sets that are involved in this work and the relationships between them are summarized in Table I.

Unless otherwise specified, matrices in the secondary AO basis are denoted by bold Roman letters (e.g., F, P), while those in the primary basis are by bold calligraphic Roman letters (e.g., ℱ, 𝒫). To concisely show the character of quantities within a nonorthogonal basis, tensorial notation will be used in the derivation, i.e., covariant (subscript) and contravariant (superscript) indices are distinguished, following Ref. 107 and the appendix of Ref. 89. For instance, a matrix element denoted by $BXαXμ$ indicates that matrix B has rows corresponding to contravariant secondary basis functions and columns corresponding to covariant primary basis functions, and these basis functions belong to the same atomic center X. Einstein summation convention is applied for contractions between contravariant and covariant indices, except for summations over different atomic centers, which will be written out explicitly.

The search for the MAB described in Sec. II B requires an inexpensively calculated reference density matrix in the secondary (target) basis. A converged SCF solution in a small PRB serves this purpose. The PRB is constructed by projecting the reference basis functions, ${|ωAα̃〉}$⁠, into the space spanned by the secondary basis on each atom:⁶⁶ 

Here, $S12AνAα=〈ωAν|ωAα̃〉$ is the overlap between the (unprojected) reference basis and the secondary basis functions, while S_A is the overlap metric of secondary basis functions on atom A. Throughout this paper, B is used to denote atom-blocked matrices containing contraction coefficients of the secondary basis functions on each atomic site, which defines a primary basis.

Since the reference basis is small (to be specified later) while the secondary basis is close to the CBS limit, the block in B_ref will be very sparse since the contraction coefficients for the high angular momentum components of the secondary basis all vanish during the projection procedure. At this stage, an SCF calculation is performed in the PRB, by solving the following generalized eigenvalue equation:

where ℱ and 𝒮 have the dimension of the PRB, and they can be transformed from their counterparts in the secondary basis using the B_ref matrix,

In reverse, the PRB density matrix, 𝒫, can be projected into the secondary basis via the following transformation:

Since the PRB is an exact subset of the secondary basis, no information in 𝒫 is lost upon projection into the latter (Eq. (4)), i.e., P_ref and 𝒫 contain the same information about the chemical environment (this is not true if 𝒫 is optimized in the reference basis directly without doing the projection). We call this special property of P_ref “PRB-representability.” The final PRB density matrix 𝒫 therefore becomes the reference used in the search for the MAB.

Sec. II A employs a basis defined by a fixed atom-blocked transformation (the PRB) and converges a density matrix in it. With the fixed P_ref in hand, our goal now is to optimize an energy-like function with respect to a variable B matrix that defines the MAB. (Note: In the following discussion B exclusively denotes the MAB coefficients.) Since a single diagonalization minimizes Tr[PF] for a chosen number of electrons⁶⁴ when F is given, we shall, by analogy, minimize $Tr[P̃F]$⁠, where $P̃$ is a “MAB-representable” density matrix in the secondary basis,

𝒟 is a density matrix in the MAB that derives from the fixed P_ref.

However, the MAB has smaller rank than the PRB, and the spaces spanned by them are rather different, so P_ref will not be exactly representable by the MAB. There exist many possible ways to construct $P̃$⁠. We choose to project P_ref into the space spanned by the MAB first, which gives 𝒟, then transform it back into the secondary basis,

where σ refers to the overlap metric of the MAB,

Recognizing that the projector into the MAB space is

our surrogate energetic objective function becomes

For brevity, in the following derivation we use P instead of P_ref for the fixed reference density matrix.

We note that the MAB-representable density matrix $P̃$ usually does not contain exactly the right electron count. While the exact N_elec is given by Tr[PS], utilizing the idempotency of R (based on Eqs. (8) and (9), it is straightforward to derive RSR = R), we have

The inequality arises because the reference density matrix P is usually not MAB-representable.

Gradient-based optimization can locate the optimal B as the minimizer of Eq. (10). The initial guess for the MAB (and its orthogonal complement in the span of the secondary basis functions on each atomic site, which is denoted by V) is obtained by diagonalizing atomic blocks of the reference density matrix, appropriately transformed⁷⁵ as $PA′=XATSAPASAXA$⁠, where S_A is the overlap matrix of the secondary basis functions on A, and X_A is the canonical orthogonalizer for them. With U_A representing the eigenvectors of $PA′$⁠, the initial B and V are set to

where $(UA)AiAp$ denotes the eigenvectors corresponding to the m_A largest eigenvalues of $PA′$ (m_A is the rank of minimal basis for atom A), and $(UA)AaAp$ denotes the remaining eigenvectors. This gives the initial partitioning of the Hilbert space that can be represented symbolically as

Since the MAB functions (and the complementary ones) are constructed by on-site contractions of the secondary basis, the variables that parameterize the MAB are intra-atomic orbital rotations. Akin to Ref. 52, a single on-block unitary transform is parameterized by the exponential of an antisymmetric matrix,¹⁰⁸ which ensures that the updated atomic orbitals stay on the same manifold,

where $CXrXμ$ denotes the union of the MAB and the complementary functions on atom X. To enforce antisymmetry of θ, it is further parameterized by Δ which contains all the independent variables,

The desired gradient, evaluated at Δ = 0, is

E is invariant with respect to orbital rotations within the MAB space (p = i, q = j), or within the space of complementary excluded vectors (p = a, q = b), as these rotations leave R unchanged. Therefore, the non-zero gradient comes only from variations of Δ^ZiZa. Using the identities

and

the desired gradient expression is given by

where, for brevity, G = ∂E/∂R as defined in Eq. (16). More details about the derivation of Eq. (19) is provided in  Appendix A.

Once the gradient at the current position is computed, the optimization algorithm will generate a new step (Δ) based on it (and the previous gradients and steps). The equations for the exponential transformation were derived in Ref. 108. The update for the MAB can be represented as

U and p stand for eigenvectors and eigenvalues of the matrix quantity ΔΔ^†, respectively, and note that the unitary transformations are atom-blocked operations. When the iterative optimization converges, B represents a minimal basis energetically adapted to the chemical environment described by the reference density matrix (from PRB-SCF). Figure 3 illustrates the MAB optimization procedure. Finally, we note that for unrestricted cases, the MABs for α and β electrons are optimized separately (they are completely decoupled), using the same objective function form.

With the MAB defined, a converged SCF solution can be obtained in this basis. The SCF energy in the fixed MAB will be an approximation to the energy evaluated by PAO-SCF, which directly minimizes the SCF energy with respect to the generators of the MAB as well as the variables defining the density matrix. These two approaches will be compared in Sec. IV.

The size of a minimal basis only depends on the principal quantum number of the atom’s valence shell, since a complete set of angular momentum functions are needed to fulfill the requirement of spatial isotropy. This definition usually works very well, but there are two types of exceptions. First, in some cases, the standard rank of the minimal basis includes redundant functions. For example, the minimal basis of lithium (n = 2) consists of 5 functions, although only two of them are required to describe its 1s²2s¹ configuration. The same applies to many other electron-deficient species like cations and radicals. The presence of redundant functions causes difficulties in converging the MAB optimization procedure. Second, in some cases, the standard rank of the minimal basis is too small to accurately describe the bonding. Examples include some hypervalent molecules (e.g., SO₃ and ClF₃), and, occasionally, molecular anions. In such cases, the standard rank of the MAB will lead to larger errors in the resulting molecular energies, which can be greatly reduced if a certain number of additional MAB functions are judiciously added to the appropriate atomic centers.

In both cases, we can adjust the rank of the MAB appropriately based on information that is already available from the initial PRB-SCF calculation. The resulting procedure, shown in Algorithm 1, can either truncate or augment the MAB dimension on each atom. The number of significant eigenvalues (N_Sig) for each atom is set to the number of eigenvalues of $PA′$ that are above a first threshold (thresh1, whose default value is 0.01). The MAB dimension will be reduced to N_Sig if that is smaller than a minimal basis (N_min). On the other hand, when N_min < N_Sig, the algorithm expands the MAB dimension by the number of eigenvalues beyond N_min that satisfy E_A(i)/E_A(N_min) > thresh2 (thresh2 has a default value of 0.02, i.e., eigenvalues that are larger than E_A(N_min)/50 will be included), which will allow a lower optimized MAB-SCF energy. The default values of these two thresholds are empirically determined based on the performance on the hypervalent molecules in the G2 set (see Sec. IV).

The objective function given by Eq. (10) can be rewritten as follows:

where $C̃o=RSCo$ represents the PRB-optimized occupied MOs after being projected into the MAB space. For stable species, the energies of occupied MOs should all be negative, and thus minimizing E corresponds to retaining as many of the bound electrons as possible.

A disaster occurs in the MAB optimization if an occupied MO has a positive energy, because minimization will result in loss of those electrons. With inexact functionals, this occasionally happens for anions. For example, the energy of the three 2p orbitals in F⁻ is 0.001 Eh with the B3LYP functional (hybrid), and 0.056 Eh with B97-D (pure). With these functionals, the resulting value of Tr[RSPS] (number of electrons captured by $P̃$⁠) is close to 2 when the MAB is optimized, which indicates that the six 2p electrons are missing! This is completely unphysical, and causes the SCF energy computed with the MAB to be qualitatively incorrect.

Such difficulties can be avoided by modifying the eigenvalue structure of F to ensure that all occupied levels are negative. This can be done by applying a uniform shift to all the eigenvalues,

using C^TSC = I. The shift, λ, is set to be

where 0 < α ≤ 1 so that the zero energy lies between the HOMO and LUMO calculated by PRB-SCF (note: such a shift is applied only when ϵ(HOMO) is detected to be positive). The “mixing” parameter α = 0.75 is empirically selected based on the performance of our method on small (monoatomic and diatomic) anions (for more details, see Sec. 1 of the supplementary material¹⁰⁹). Replacing F with F′ in Eq. (10) gives a modified objective function for the MAB,

When λ > 0, the new term resembles a penalty for losing electrons, which can be made explicit by adding an additional constant, λN_elec, to the RHS,

Based on the data presented in Refs. 89 and 90, a significant difference exists between PAO-SCF and exact SCF results. To reduce this gap, computationally inexpensive correction schemes based on perturbation theory are useful. Analogous to the dual-basis methods, the converged MAB-SCF solution serves as the primary basis reference, and the contribution of non-Brillouin singles to the second-order perturbative (PT2) energy correction is given by⁶³ 

Here F denotes the Fock matrix built upon the MAB-SCF density projected into the secondary basis: $F=F(P̃)$⁠. The first-order T-amplitude satisfies the following linear equation:

In the pseudo-canonicalized MO basis (obtained by diagonalizing F_OO and F_VV separately, see  Appendix C), Eq. (27) reduces to a simpler form,

Correspondingly, the perturbative energy lowering becomes

which can be interpreted as an energy-weighted steepest descent (or an approximate Newton) step.^63,90 Alternatively, other corrections that involve a full diagonalization of the Fock matrix can be applied, such as the aforementioned DB-SCF (only slightly different from PT2) and DFPC methods. The latter performs a single update of the density matrix in the secondary basis (by diagonalizing F), and then recomputes the full SCF energy based upon that density matrix (the result will thus be variational).

A pilot implementation of this new SCF scheme is accomplished in a development version of the Q-Chem 4.3 package.¹¹⁰ A preconditioned limited-memory BFGS (L-BFGS) algorithm^111,112 is implemented for solving the MAB optimization problem efficiently. The inverted on-diagonal blocks of the Hessian matrix for the objective function (second derivatives of the function value with regard to orbital rotations on the same atomic site) are employed as the preconditioner of the L-BFGS algorithm. In most scenarios, this preconditioning strategy leads to convergence of the MAB optimization in a reasonable number of iterations (10¹–10²), with only moderate additional cost to evaluate the preconditioner. More details about the preconditioned L-BFGS algorithm and the evaluation of the on-diagonal blocks of the Hessian are provided in  Appendix B.

In the current implementation, all the density matrix updates are computed by diagonalizing ℱ (Fock matrix in the dimension of PRB or MAB), and the only diagonalization in the secondary basis dimension (N_v × N_v to be exact) is required by the perturbation correction. However, to obtain ℱ (for the time being) still requires contracting the ERI tensor with a PRB- or MAB-representable density matrix in the secondary basis (using restricted KS Fock matrix as an example),

where κ is the proportion of exact exchange in the employed functional. Then F is transformed back into the primary basis through Eq. (3). This choice is actually less efficient, because quantities in the primary basis, like $Pαβ$⁠, $Cαi$⁠, can be directly utilized to construct ℱ, which will significantly reduce the dimension of the contraction. Therefore, with our preliminary implementation, we will focus on validating the accuracy of MAB-SCF (PC) in this work, and the potential efficient implementation of this method will be briefly discussed in Sec. VI.

All the results for TC and NC datasets are generated with the KS-DFT routines in Q-Chem 4.3 as well. A (75, 302) grid (75 radial shells with 302 Lebedev points in each) is used for all employed exchange-correlation (XC) functionals, and the SG-1 grid¹¹³ is used for the VV10¹⁸ NLC functional. Unless otherwise noted, 6-31+G(d) is used as the reference basis in PRB-SCF. Smaller reference bases could alternatively be employed, at the cost of diminished accuracy (see the results provided in Sec. 2 of the supplementary material¹⁰⁹). The optimization of the MAB converges to 10⁻⁶ a.u., while all the SCF calculations converge to 10⁻⁸ a.u. (RMS of the gradient). To determine the appropriate dimension of the MAB, the default values of thresh1 and thresh2 adopted by Algorithm 1 are set to 0.01 and 0.02, respectively. We note that, in this work, the “adding vector” strategy is by default turned off and only utilized for specified hypervalent molecules.

We start investigating the accuracy of our method by performing a series of preliminary tests on the G2 thermochemistry set.³⁵ To test the quality of the optimized MAB, the performance of MAB-SCF on the G2 set is compared to PAO-SCF, since the latter gives the limiting behavior of an atom-centered minimal basis. For the reasons discussed in Sec. II C, the minimal basis models (including MAB and PAO) are not sufficient for describing hypervalent molecules. Thus, we designate molecules containing Al, Si, P, S, Cl centers that are coordinated by highly electronegative atoms (e.g., O, F, Cl) as hypervalent, including SO, ClO, SO₂, AlF₃, AlCl₃, SiF₄, SiCl₄, PF₃, ClF₃, and (CH₃)₂SO, and exclude them from the test set preliminarily.

The results for this “pruned” G2 set (138 molecules) computed with three functionals are collected in Table II (aQZ is the target secondary basis). To make it a fair comparison, molecules that fail to converge their PAOs (listed in the table footnotes) are also excluded when reporting the statistical errors. At the SCF level (i.e., without PT2 correction), PAO-SCF and MAB-SCF show similar accuracy with respect to the exact SCF results for all three tested functionals. Surprisingly, the MAB-SCF results are slightly better as a result of error cancellation (PAO-SCF is exact for atomic energies). Applying the PT2 correction significantly reduces the errors of both schemes. For the two pure functionals (B97-D and B97M-V), the RMSDs of MAB+PT2 are smaller than 0.1 kcal/mol (∼0.05 kcal/mol), and they are close to those of PAO+PT2. For B3LYP, the RMSDs of both schemes noticeably increase, which suggests diminished effectiveness of PT2 when hybrid functionals are used. Nevertheless, we notice that the performance of MAB+PT2 is rather similar to that of PAO+PT2.

To better compute the energies of the hypervalent molecules, we increase the dimensions of their MAB based on Algorithm 1. The modified MAB+PT2 results are compared with those using the standard minimal basis dimensions in Table III. For a majority of these molecules (ClO, SiF₄, SiCl₄, PF₃, (CH₃)₂SO, and presumably SO₂ and ClF₃), the errors are reduced by over 10 times by using the “adding vector” strategy. The degree of inadequacy of the conventional minimal basis dimension is perhaps a measure of molecular hypervalency. Indeed, species like SiF₄, PF₃ do not formally violate the “octet” rule, which indicates that molecular hypervalency may exist beyond its usual definition. On the other hand, AlCl₃ and AlF₃ do not seem to be typical hypervalent species, because the use of standard minimal basis dimensions does not result in errors as large as those for the other molecules listed in Table III.

Combining these specially treated hypervalent molecules with the other 138 molecules computed with the standard MAB model, the overall RMSD of MAB+PT2 against conventional SCF results for the G2 set is 0.033 kcal/mol. This result is only minimally different from the RMSD for the “pruned” G2 set (with B97M-V), by contrast with the poor results for standard MAB+PT2 and PAO+PT2 when hypervalent molecules are included. Therefore in the later G2 tests, unless otherwise specified, statistical errors evaluated including with all the molecules will be reported with the hypervalent ones separately treated via Algorithm 1.

Table IV shows how the dimensions of the MAB are increased on the central atoms of the hypervalent molecules after applying Algorithm 1. According to the rightmost column, the modified MAB function counts are usually close to those of a minimal basis plus one set of d (polarization) functions. Although AlF₃ and AlCl₃ do not show strong hypervalent character according to Table III, the Al atom nonetheless gains additional MAB basis functions.

In general, the MAB optimization problem is considerably easier to converge than the aforementioned “double” PAO-SCF optimization. In contrast to PAO-SCF that encounters several convergence problems (mentioned in Table II), no MAB convergence failure is detected for the entire G2 set with all three tested functionals. Furthermore, in contrast to PAO-SCF, the MAB optimization is decoupled from density matrix optimization, and thus MAB iteration counts will not directly affect the required number of SCF cycles. This is extremely important because of the much more significant cost per iteration for the latter. Table V lists molecules in the G2 set that require over 100 iterations to converge their PAOs (using B97M-V/aQZ). Apart from several convergence failures, the PAO-SCF scheme requires an enormous number of Fock matrix constructions for some molecules, such as NaCl. The MAB scheme, on the other hand, attains the optimized adaptive basis in fewer iterations for most of these molecules. Even in cases like COS and C₂H₄O (oxirane) where the iteration counts for optimizing the adaptive basis are similar, MAB-SCF is still far more efficient because many fewer Fock builds are required. Therefore, MAB-SCF appears to be a more feasible adaptive basis SCF scheme than PAO-SCF, with comparable accuracy.

Table II already suggests that the performance of MAB+PT2 is not completely functional-independent. A clear difference exists between using pure and hybrid functionals. Therefore, we must investigate the performance of this method when different flavors of density functionals are employed. Using the G2 set with aug-cc-pVQZ as the secondary basis, Table VI explores the performance of 13 density functionals. The first seven functionals do not contain exact exchange, including three GGAs (B97-D, BLYP,^26,27 PBE⁵) and four meta-GGAs (TPSS,⁸ MGGA_MS1,¹¹⁵ M06-L,¹¹⁶ B97M-V). For these functionals, the MAB+PT2 scheme demonstrates good accuracy (M06-L is the largest outlier), while the RMSDs computed by MAB+DFPC are roughly twice as large. Thus PT2 appears preferable.

The other six functionals in Table VI are hybrid functionals, including TPSSh (10%),¹¹⁷ B3LYP (20%), PBE0 (25%),¹¹⁸ M06-2X (54%), ωB97X-D (RSH),¹¹⁹ and ωB97X-V (RSH) (“%” denotes the proportion of exact exchange). They show results that contrast with those of the pure functionals. For these hybrids, PT2 undershoots the exact SCF energy, and the size of the errors grows roughly with the amount of exact exchange. This leads to unsatisfactory accuracy (RMSD > 0.2 kcal/mol) for functionals that contain more exact exchange than B3LYP. The MAB+DFPC approach, on the other hand, shows comparatively better performance across the hybrids (RMSDs are around 0.1 kcal/mol), than PT2. For the two RSH functionals, the results of MAB+DFPC are 5–6 times more accurate than MAB+PT2.

We conclude that DFPC should be used as the correction to MAB-SCF when hybrid functionals are employed, at the expense of one more Fock build in the secondary basis. On the other hand, MAB+PT2 is less expensive and more accurate for pure functionals.

Different basis sets are employed for KS-DFT calculations, based on considerations such as accuracy, efficiency, and user experience. With the B97M-V functional, we assess the performance of MAB+PT2 with several widely used basis sets, including those in Dunning’s correlation-consistent series (aTZ, QZ, aQZ),^28,29 Jensen’s polarization-consistent series (apc-2, pc-3, apc-3),^120–122 and the Karlsruhe “def2” series (TZVPPD, QZVPP, QZVPPD).³³ The popular “large Pople” basis set 6-311++G(3df,3pd)^123,124 is also included. These basis sets are augmented triple-ζ quality or higher, because our goal is to approach CBS limit results.

The RMSDs for the G2 set using different target secondary basis sets are displayed in Figure 4. The errors are typically below 0.1 kcal/mol, which indicates excellent transferability. The best performance is achieved by two quadruple-ζ basis sets with diffuse functions, aQZ and QZVPPD. For their unaugmented counterparts, QZ and QZVPP, the RMSDs are slightly larger, although still satisfactory. Since diffuse functions typically have no major impact on the accuracy of evaluated energetics for bonded interactions, the compatibility of MAB+PT2 with these unaugmented basis sets is helpful. The largest RMS errors are produced by Jensen’s pc-3 and aug-pc-3 basis sets. This is mostly due to the inaccurate atomic energy for Li: the RMSDs are reduced to 0.063 (pc-3) and 0.059 (apc-3) kcal/mol if we exclude three Li-containing molecules (Li₂, LiF, and LiH) from the G2 set. As these are high-quality basis sets, the outlier might be due to poor compatibility of pc-3/apc-3 with the employed reference basis, 6-31+G(d), for the Li atom. Computational cost aside, any of the largest basis sets approach the true CBS limit and also the top accuracy of our method.

Ultimately, MAB-SCF with a perturbation correction scheme (MAB-SCF (PC)) might replace conventional SCF in KS-DFT calculations to approach the CBS limit. Therefore, we want to see how MAB-SCF (PC) approaches the CBS limit with the increasing size of the secondary basis. For that purpose, we extend Figure 5 (which motivates this work) with the functional RMSDs evaluated by MAB+PT2 (for B97-D, B97M-V) and MAB+DFPC (for B3LYP, M06-2X, ωB97X-V), as shown in Figure 5. These results are quite encouraging: the conventional SCF convergence towards the CBS limit of each functional is closely reproduced by MAB-SCF (PC). At the aQZ level, the largest difference between MAB-SCF (PC) and exact SCF results is only about 0.03 kcal/mol (for B97-D and B3LYP), which is below 1% of the intrinsic (CBS) error of the functional itself, and thus negligible. Larger differences between these two sets of results exist at the aDZ level (the largest gap is 0.26 kcal/mol for ωB97X-V/aDZ), but this is not relevant to our target of the CBS limit.

We will assess how the performance of MAB-SCF (PC) transfers to other thermochemistry (TC) datasets. Three density functionals (B97-D, B97M-V, B3LYP) will be employed to examine the accuracy of MAB-SCF (PC), using aQZ as the secondary basis, and the perturbation correction schemes are applied in the same way as in Figure 5. First we consider the W4-11 dataset,¹²⁵ which includes 99 bond dissociation energies (BDE99), 707 heavy-atom transfer reaction energies (HAT707), 20 isomerization energies (ISO20), 13 nucleophilic substitution reaction energies (SN13), and 140 total atomization energies (TAE140). Note that the multi-reference (MR) species in W4-11 are included in this test. Thirteen species are separately treated as hypervalent molecules: AlF₃, AlCl₃, SiF₄, P₄, SO, SO₂, SO₃, S₂O, S₂, S₃, S₄, ClO, and OClO (see the previous discussion for the G2 set).

Table VII contains the RMSDs of the MAB-SCF (PC) approach (against exact SCF) for W4-11, where different TC categories have been separated. The overall performance is similar to that for the G2 set, which shows encouraging transferability. Taking B97M-V as an example, the smallest and largest RMSDs of MAB+PT2 (vs. exact SCF) are obtained on SN13 (0.02 kcal/mol) and HAT707 (0.06 kcal/mol), respectively, while the corresponding functional RMSDs (vs. W4 reference) are 1.39 kcal/mol and 3.90 kcal/mol. Therefore, the errors caused by replacing conventional SCF with MAB-SCF (PC) are usually smaller than intrinsic functional errors by one or two orders of magnitude. A more straightforward comparison is provided by Figure 6: for B97M-V, the functional RMSDs computed via the MAB+PT2 approach show almost no difference compared to those by normal SCF method.

We also examined the performance of MAB-SCF (PC) on other TC datasets, including adiabatic ionization potentials and electron affinities (21 for each: G21IP and G21EA),^126,127 38 non-hydrogen transfer and 38 hydrogen transfer barrier heights (NHTBH38¹²⁸ and HTBH38¹²⁹), and 14 alkane isomerization energies (Pentane14¹³⁰). The computational details are identical to those for W4-11, except for the anions in G21EA, HTBH38 and NHTBH38, where the modified MAB objective function (introduced in Sec. II D) is automatically applied to avoid the emergence of unphysical results. The RMSDs of MAB-SCF (PC) are presented in Table VIII. We see that HTBH38 and Pentane14 are relatively easier cases for MAB-SCF (PC) to approximate the exact SCF result; and for G21IP and NHTBH38, the size of the RMSDs is similar to that for those W4-11 subsets (e.g., BDE99 and SN13).

The largest RMSD appears on the G21EA dataset. Although none of the results are qualitatively incorrect by applying the modified MAB objective function, there are several molecular anions whose absolute energies evaluated by MAB-SCF (PC) are rather unsatisfactory: NO⁻, PO⁻, $O2−$⁠, and $S2−$⁠. Based on the discussions in Sec. II C, the “adding vector” strategy may also be applied for these electron-abundant species. As a result, the absolute energies of these anions are significantly improved, as shown in Table IX. With these 4 molecular anions specially treated, the RMSDs for the G21EA dataset are recalculated and the results are also presented in Table VIII, which turn out to be more comparable to the RMSDs for other TC test sets.

One of the key improvements in modern density functionals is for non-covalent interactions (NC). Therefore, we assess the performance of MAB-SCF (PC) on several NC datasets, including A24 (24 small NC complexes),³⁰ S22 (22 diverse small- to medium-sized NC complexes at the equilibrium geometries),^131,132 HB15 (15 ionic hydrogen bond interactions),¹³³ H2O6Bind8 (binding energies of eight configurations of water hexamers),^134,135 and FmH2O10 (binding energies of 10 configurations of F⁻(H₂O)₁₀^134,135). Since many NC interactions have smaller magnitudes than TC energy differences, and some modern density functionals are able to achieve very small errors for them (e.g., B97M-V’s unweighted RMSD for 1458 non-covalent interactions is 0.22 kcal/mol²⁰), higher accuracy is needed for the MAB approach to match the exact SCF results.

As before, three density functionals (B97-D, B97M-V, B3LYP-D3(0)) are employed to assess the performance of MAB-SCF (PC) on these NC datasets. To avoid counterpoise (CP) corrections, we choose def2-QZVPPD as the secondary basis, which has fewer functions than aQZ but generates even smaller BSSEs.¹³⁶ Table X contains the resulting RMSDs (vs. exact SCF) for these NC datasets. Very small MAB-SCF (PC) errors are found for dimer binding energies (A24, S22, and HB15), with all three functionals. This is very encouraging for treating the most common non-covalent interactions by MAB-SCF (PC) instead of conventional SCF. Larger errors appear for the cluster binding energies (H2O6Bind8, FmH2O10), due to the larger magnitude of these interactions (H2O6Bind8: −40 to −50 kcal/mol; FmH2O10: about −170 kcal/mol) and the uniformity of their interaction types (systematic errors accumulate on one single direction through all the data points). The largest outlier appears on FmH2O10 when the B3LYP-D3(0) functional is employed, mostly due to the less accurate MAB-DFPC result for F⁻ (the error is 0.3 kcal/mol vs. exact SCF) with that functional.

The accurate description of these non-covalent interactions by B97M-V can also be reproduced using MAB-SCF (PC), and the resulting RMS errors (vs. reference data) are compared with the exact SCF results in Figure 7. Only minimal differences exist between two sets of RMSDs for the dimer binding energies, while the monotonic deviations of the MAB-SCF (PC) results for each single data point contribute to more pronounced differences for the clusters. Nevertheless, even for H2O6Bind8 where the largest deviation occurs, the RMSD vs. exact SCF results (0.085 kcal/mol) is only 0.1%–0.2% of the magnitude of the corresponding binding energies. In terms of the evaluation of relative energies, fairly insignificant errors are caused by using MAB-SCF (PC).

Finally we test the accuracy of MAB-SCF (PC) on large non-covalent complexes, using the seven dispersion-bound systems of the recently proposed L7 dataset,¹³⁷ including (the abbreviations simply follow Ref. 137): stacked circumcoronene-adenine dimer (C3A), stacked circumcoronene with a Watson-Crick G-C base pair (C3GC), parallel displaced coronene dimer (C2C2PD), stacked Watson-Crick G-C base pairs (GCGC), stacked guanine trimer (GGG), parallel stacked octadecane dimer (CBH), and phenylalanine residue trimer (PHE). Due to the tremendous computational effort required for these systems, our calculations were performed with two pure functionals (B97-D, B97M-V) using aug-cc-pVTZ as the secondary basis. Binding energies of these complexes evaluated by conventional SCF and MAB-SCF (PC) are compared in Table XI. For most of them, the differences between the results given by two SCF schemes are much smaller than 1% of the magnitude of their binding energies, which indicates the satisfactory accuracy of MAB-SCF (PC) for these large-scale non-covalent interacting systems. The only exception is the guanine trimer (GGG), which is primarily due to its very weak binding energy (−2.33 kcal/mol by B97M-V/aTZ). Due to the lack of reliable reference values for the time being, we cannot meaningfully assess the intrinsic functional errors for L7.

The tests performed on a broad range of systems demonstrate that the overall accuracy of this new approximate SCF method seems to be encouraging. As the CBS limit is approached, the functional RMSDs for several representative TC and NC datasets (G2, W4-11, S22) evaluated by MAB-SCF (PC) show very small differences compared to the exact SCF results, as demonstrated in Figures 5, 6, and 7. In particular, “B97M-V/QZVPPD/MAB+PT2” (or “B97M-V/aQZ/MAB+PT2”) turns out to be a promising model chemistry, because of B97M-V’s very good accuracy for both TC and NC, its moderate computational cost as a semi-local functional, and its compatibility with the MAB+PT2 approach.

The MAB-SCF (PC) procedure is similar to that of dual-basis methods (including DB-SCF and DFPC), since they both first compute an SCF solution in a primary basis (preferably a subset of the secondary basis) then correct it using perturbation theory. It is encouraging that MAB-SCF (PC) is able to achieve comparable accuracy to DB-SCF with a primary basis adaptively prepared on the fly, whose size is also much smaller than those well-trained DB-SCF basis subsets⁶⁷ when the secondary basis approaches the CBS limit. Some comparisons of these two approaches and related discussion are provided in Sec. 4 of the supplementary material.¹⁰⁹ 

At least three challenges remain in successfully applying the MAB-SCF (PC) method. First is the challenge of hybrid density functionals. According to Table II, there is no degradation in terms of the quality of the MAB when a hybrid functional is used: B3LYP in fact has the smallest RMSD for MAB-SCF among the three tested functionals. The reason for the poor performance of MAB+PT2 for hybrid functionals should therefore reside in the PT2 correction represented by Eq. (29). Assuming that the MABs calculated by two functionals are of similar quality, the numerators in Eq. (29) should also have fairly similar sizes. Thus, functional-dependent differences in the magnitude of the PT2 correction will be largely determined by the denominators (orbital energy differences) in Eq. (29). It turns out that the smaller gaps calculated by pure functionals help reduce the errors of MAB-SCF more effectively, while hybrid functionals that usually give larger orbital energy gaps undercorrect. This numerical result is intriguing since the gaps computed by the hybrids are usually deemed to be physically more correct. Although the accuracy of MAB+DFPC is satisfactory in most tested cases, an extra Fock build using the density matrix in the secondary basis (not MAB-representable) is required, which would be preferable to avoid (the reason will be elucidated below). Additionally, based on Table VI, the accuracy of MAB+DFPC slightly degrades for RSH functionals as well despite its clear advantage over MAB+PT2.

A second limitation is that hypervalent molecules are treated separately in this work by modifying the MAB dimensions in a semi-automated way (choosing to use the “adding vector” approach is user-specified). This is because our algorithm and its related thresholds were chosen to attain the desired accuracy for hypervalent molecules. This option increases the size of the MAB (and thus the computational cost of MAB-SCF) and potentially degrades the convergence of MAB optimization (see the discussion in Sec. 3 of the supplementary material¹⁰⁹). Further work on wisely adjusting the algorithm parameters, or possibly refining the present scheme, is desirable.

The third general challenge is efficient computational implementation. Although this paper is mainly about the formulation MAB-SCF (PC) and the validation of its accuracy, the ultimate goal is to serve as an inexpensive substitute for conventional SCF as the CBS limit is approached. Compared to normal SCF, our method involves far fewer degrees of freedom. The cost of a single density matrix update, if achieved by diagonalizing the Fock matrix, can be reduced by up to a factor of (n/m)³, where n is the size of the secondary basis, and m is the size of the PRB or MAB. For example, given a model system (CH₂)_n, there will be 129 basis functions per –CH₂– unit if def2-QZVPPD is employed as the secondary basis, which is roughly 6 times as large as the PRB (22 functions per unit), and about 18 times as the MAB (7 functions/unit). Therefore, the prefactor of the computational cost of this cubic scaling step is significantly reduced in our scheme. A linear equation solve in the secondary basis dimension is still required to obtain the PT2 correction, while that is far less expensive than diagonalization.

Alternatively, $O(N)$ scaling electronic structure methods can be potentially applied to the MAB-SCF step, since the overlap matrix of the MAB is extremely well-conditioned, as demonstrated by Table XII. The diagonalization-free density matrix update algorithms introduced by Ref. 52 could be used for the MAB-SCF. However, the feasibility of these methods for PRB-SCF is not so clear because one set of diffuse functions is contained in 6-31+G(d).

In practice, however, even for systems as large as the L7 complexes, the Fock matrix construction step still dominates the computational cost of each SCF cycle due to its large prefactor when high accuracy is sought. This is despite the fact that asymptotically that step scales quadratically¹³⁸ (or even linearly with special algorithms)^40–47 with respect to system size. In our pilot implementation, the Fock matrix is built from PRB- or MAB-representable density matrices that are of secondary basis dimensions, which costs almost the same as a Fock build in a conventional SCF calculation within the secondary basis. Therefore, to speed up these Fock build steps by taking advantage of the properties of the PRB and MAB becomes the most urgent task for MAB-SCF (PC) to outperform conventional SCF in terms of computational efficiency, especially for medium-sized systems. PRB-SCF can be reformulated as a conventional SCF calculation within a basis set whose size and shell structure are identical to 6-31+G(d), due to the elimination of high angular momentum functions during the basis set projection procedure, and the cost of the involved Fock builds thereby can be significantly reduced without requiring more sophisticated techniques. For MAB-SCF, instead of using Eq. (30), the MAB Fock matrix can be constructed by forming the ERI tensor in the primary basis first then contracting it with the MAB density matrix, as

Furthermore, quantities like (αβ|γδ) can be efficiently computed (and stored) using the RI approximation, due to the highly compact and atom-blocked structure of the MAB. We intend to present an optimized implementation of MAB-SCF (PC) and the resulting timings vs. conventional SCF as soon as these challenges are adequately addressed.

In this work, we proposed a new minimal adaptive basis (MAB) SCF method that can be used in KS-DFT calculations. Its objective is to permit approach to the complete basis set (CBS) limit without explicitly performing the calculation in very large AO basis sets. The key aspects of this paper can be summarized as follows:

• An MAB is obtained as a molecule-adapted, atom-blocked transformation from the secondary AO basis. We have developed a viable optimization method to obtain an MAB using an inexpensive energy-like surrogate function based on a reference SCF calculation in a projected basis of moderate size (the PRB). Compared to exact energy optimization with respect to the transformation (the PAO method), our MAB yields similar total energies with far fewer convergence issues.

• A preconditioned L-BFGS algorithm that requires the gradient and the on-diagonal blocks of the Hessian of the objective function is implemented to solve the MAB optimization problem. In addition, an approach that modifies the MAB dimension based on chemical environment is proposed, and demonstrated to be necessary for hypervalent molecules. These ideas can potentially be used in PAO calculations and other adaptive basis models as well.

• Perturbation corrections (PC) are applied to MAB-SCF to approach the desired accuracy. This resembles DB-SCF without the need to select or develop the paired basis subset. The resulting accuracy is assessed on numerous TC and NC datasets. Measured against exact SCF results, MAB-SCF (PC) generates <0.15 kcal/mol RMSDs for most of the tested TC datasets, and even smaller errors (usually <0.1 kcal/mol) for the NCs. Encouragingly, as the CBS limit is approached, the MAB-SCF (PC) method deviates from full SCF one or two orders of magnitude less than the inherent errors in today’s best functionals.

• Future work includes further refining the MAB-SCF (PC) model, and developing an efficient implementation, as discussed in Sec. VI. We note that hybrid functionals require a different PC to achieve adequate accuracy, and that the strategy to treat hypervalent species is not yet fully automated. Nonetheless, based on the accuracy demonstrated here, and the potential computational advantages of using an atom-centered minimal basis, we believe that this approach merits further development. We hope to report further progress in due course.

This work was supported by Grant No. CHE-1363320 from the U.S. National Science Foundation. C.K.S. would like to acknowledge support from EPSRC Grant No. EP/K039156/1.

Using the same notations as in Sec. II B, the most general form of the gradient of the MAB objective function is given by

where G = SPSRF + FRSPS. With the parameterization of B by Equations (14) and (15), we have

and correspondingly the gradient becomes

When Δ stands for the orbital rotations within the MAB space, i.e., p = i, q = j, matrix C reduces to B, and the gradient vanishes because (I − RS)B = 0. On the other hand, orbital rotations within the complementary space of the MAB (p = a, q = b) have no effect on the objective function value either, simply due to the enforced on-atom orthogonality (σ_ZaZj = 0). Therefore, the only non-zero block of this gradient is resulted from the rotations between these two subspaces. If we set p = i, q = a, only the first term in Eq. (A3) remains based on the arguments above, which immediately leads to the gradient represented by Eq. (19).

The basic idea of L-BFGS is to construct the approximate Hessian (inverse Hessian in practice) for the current iteration with gradients and displacements computed in the most recent m steps, where m is the user-specified subspace size (number of “memorized” steps). If we define the gradient and displacement at kth iteration as g_k and s_k, and y_k = g_k+1 − g_k, the kth approximate inverse Hessian can be evaluated as¹¹¹ 

where

In practice, a “two-loop” algorithm which only requires evaluating vector-vector products is implemented to compute H_k acting on g_k. $Hk0$ is the preconditioner for the approximate inverse Hessian. By default, a constant scaling factor is used for $Hk0$⁠, which is considered as the unpreconditioned case here. Once $Hk0$ contains more information about the true inverse Hessian, the step generated by L-BFGS becomes closer to a Newton step, which can presumably accelerate the convergence.

If we still denote G = SPSRF + FRSPS, the Hessian of the MAB objective function (Eq. (10)) can be formally represented as follows (only on-block mixings are allowed):

where R^Δ and R^ΔΔ stands for first- and second-order derivatives of R with respect to Δ. The first four terms on the RHS of Eq. (B3) come from the “R^ΔR^Δ” term, while the rest from the “R^ΔΔ” term. The explicit forms of the involved matrix elements are

and

In practice, we also have σ_XiXj = δ_ij since on-block orthogonality is enforced. More details about the Hessian derivation can be found in Ref. 139, which carefully derived the orbital Hessian for SCF-MI (a rather similar optimi zation problem).

The preconditioner we apply to the L-BFGS algorithm is the inverted on-diagonal blocks of the Hessian, i.e., the inverse of H_XX for all the different atom blocks (X). Within the “two-loop” implementation, $Hk0$ acts on vector v = V_k−mV_k−m+1⋯V_k−1g_k, which can be divided into contributions from each atom block. Therefore, the application of the preconditioner is equivalent to solving the following linear equation on each atom block:

where u is the preconditioned vector: $u=Hk0v$⁠. Based on the property of the Hessian matrix (symmetric positive-definite), a preconditioned conjugate-gradient (CG) algorithm is implemented to solve Eq. (B8) iteratively on each atom block. The implemented preconditioner for CG is actually the inverse of the on-diagonal part (X = Y) of the last two terms in Eq. (B3).

Once MAB-SCF converges, a Fock matrix in the secondary basis can be built upon the MAB density matrix projected into the secondary basis,