Automatic Generation of Auxiliary Basis Sets - ACS Publications

Article pubs.acs.org/JCTC

Automatic Generation of Auxiliary Basis Sets Georgi L. Stoychev, Alexander A. Auer, and Frank Neese* Max Planck Institute for Chemical Energy Conversion, Mülheim an der Ruhr 45470, Germany ABSTRACT: A procedure was developed to automatically generate auxiliary basis sets (ABSs) for use with the resolution of the identity (RI) approximation, starting from a given orbital basis set (OBS). The goal is to provide an accurate and universal solution for cases where no optimized ABSs are available. In this context, “universal” is understood as the ability of the ABS to be used for Coulomb, exchange, and correlation energy fitting. The generation scheme (denoted AutoAux) works by spanning the product space of the OBS using an even-tempered expansion for each atom in the system. The performance of AutoAux in conjunction with different OBSs [def2-SVP, def2-TZVP, def2-QZVPP, and cc-pwCVnZ (n = D, T, Q, 5)] has been evaluated for elements from H to Rn and compared to existing predefined ABSs. Due to the requirements of simplicity and universality, the generated ABSs are larger than the optimized ones but lead to similar errors in MP2 total energies (on the order of 10−5 to 10−4 Eh/atom).

1. INTRODUCTION The evaluation of two-electron repulsion integrals (ERIs) is a time-consuming step in all ab initio electronic structure calculations and is often a significant bottleneck for the modeling of large systems due to its formal O(N4) scaling with the total number of basis functions (note that in practice, the scaling is better due to the use of integral screening procedures1−5). For this reason, methods have been developed for the approximate treatment of ERIs, which aim to achieve sufficient accuracy at a reduced cost. The most common of these is the resolution of the identity (RI) approximation in which products of the orbital basis functions (OBFs) {χμ(r)} (typically atom-centered Gaussian functions) are expanded over a set of auxiliary basis functions (ABFs) {φP(r)}:6 χμ (r)χν (r) ≈

∑ CμνP φP(r)

The obvious drawback of all RI methods is the requirement of an auxiliary basis set (ABS), which must provide a good fit for the product basis |μν) and therefore in principle depends on the orbital basis set (OBS). Naturally, the larger the ABS, the less significant the savings in computational time. There are several approaches to constructing the ABS. One of them involves a Cholesky decomposition (CD) of the ERIs, resulting in an ABS consisting of products of OBFs with a single parameter, the decomposition threshold δ, which controls the size and accuracy of the ABS and removes redundant functions.9 A variant, denoted aCD,10 whereby the ABFs are restricted to products of OBFs centered on the same atom, has also been proposed and is similar to the linear combination of atomic distributions (LCAD) method by Ten-no and Iwata.11,12 The altogether more popular approach is to use ABSs that have been previously optimized to be compact while producing small errors for a certain property, usually the total energy or the root-mean-square deviation (RMSD) in the ERIs themselves. It is important to note that these ABSs are designed to fit only certain contributions to the energy, e.g., those occurring in the Coulomb (usually denoted RI-J) or the Coulomb and exchange (RI-JK) contributions in HF or DFT or two-electron integrals in post-HF methods (RI-C). This approach requires fewer ABFs in a given set but necessitates the use of different ABSs for various parts of a calculation, e.g., RI-JK for HF and RI-C for post-HF methods. It is also possible to approximate the Coulomb energy using RI and employ a different approach for the exchange energy, as done in the RIJCOSX scheme.13

(1)

P

CPμν

where the expansion coefficients are best obtained in the Coulomb metric,7,8 i.e., by minimizing the self-repulsion of the residual (Rμν|Rμν), where Rμν(r) = χμ (r)χν (r) −

∑ CμνP φP(r) P

(2)

Thus, the ERIs (in Mulliken notation) are approximated as (μν|κλ) ≈ (μν|κλ)RI =

∑ (μν|P)(P|Q )−1(Q |κλ) PQ

(3)

where μ, ν, κ, and λ label the OBFs, and P and Q are the ABFs. This approximation allows the formal scaling of the calculation of Coulomb integrals in DFT to be reduced to O(N3), while more modest but significant savings in computational time are achieved for HF and post-HF methods, where the exchange contribution and ERIs involving virtual orbitals are also approximated. © 2016 American Chemical Society

Received: October 25, 2016 Published: December 22, 2016 554

DOI: 10.1021/acs.jctc.6b01041 J. Chem. Theory Comput. 2017, 13, 554−562

Article

Journal of Chemical Theory and Computation

Another example are basis sets tuned for specific properties such as those developed by Rappoport and Furche.37 The authors used the unaltered ABSs with their augmented OBSs and found the RI errors to be negligible in their tests, but that is not necessarily true in the general case. Furthermore, it is often required to assess the accuracy of the RI approximation, ideally by comparison against a calculation without RI. However, when the latter is prohibitively expensive, a larger ABS may be used. Therefore, it is useful to have a robust and universal solution which guarantees accuracy yet still exploits the efficiency of RI. This was the motivation behind the development of the AutoAux procedure described in this paper. In the following section, we give the principal ideas behind the generation scheme and set out a detailed algorithm of the procedure. Next, a series of benchmark results are put forward to assess the accuracy and efficiency of AutoAux followed by final remarks and an outline of possible future improvements.

Procedures for the optimization of ABSs have been developed and employed in the Karlsruhe group to produce a number of RI-J,14−16 RI-JK,16−18 and RI-C19−22 ABSs for different parent OBSs and for elements from H to Rn (excluding RI-C ABSs for the Lanthanides). A brief (and simplified) summary of these procedures is given below; for a detailed description, the reader is referred to the original publications. 1. First, a construction scheme is chosen, which defines the maximum angular momentum in the ABS, the number of functions Nl of a given angular momentum l, and a set of parameters to be optimized. These can directly be the exponents {αi} of the ABFs or the parameters αl1, βl, and γl used in an expansion formula:14,15,18 ⎛ ⎛ i − 2 ⎞2 ⎞ αil = β l αil− 1⎜⎜1 + γ l ⎜ l ⎟ ⎟⎟ , i = 2, ..., N l ⎝ N + 1⎠ ⎠ ⎝

(4)

2. AIM, CONCEPT, AND ALGORITHM OF THE AUTOAUX GENERATION PROCEDURE It is our goal to provide a scheme that allows a suitable ABS to be generated starting only from a given OBS. Our priorities are, in order of importance, as follows: 1. Accuracy. Errors in energy must be on par with those produced by optimized ABSs. 2. Universality. The ABS has to be applicable to all RIaccelerated methods and must therefore accurately fit all types of ERIs (those occurring in DFT, HF, and post-HF methods). 3. Simplicity. The ABS should be usable without optimization of exponents and contraction coefficients. For this reason, we chose an uncontracted eventempered design. 4. Efficiency. Due to the above requirements, we cannot expect the generated ABSs to be as efficient as the optimized ones, but they must be small enough to allow for a significant speedup due to the RI approximation. It should be noted that the requirements for JK-fitting and Cfitting basis sets are quite different.23 Coulomb and exchange integrals in HF and DFT contain only occupied MOs and therefore the errors in energy are fairly independent of the OBS, as demonstrated for the “universal” fitting basis sets by Weigend (referred to here as def2-J and def2-JK).15,17 For the same reason, functions with high angular momenta are not required in a JK-fitting basis set, whereas they are necessary in C-fitting basis sets because correlation also involves virtual orbitals featuring high angular momentum functions that describe angular correlation effects. Additionally, correlated calculations often employ a frozen-core approximation, in which case the ABS does not need to accommodate the core MOs. Therefore, it is computationally more efficient to devise different procedures to generate J- (for pure DFT), JK- (for hybrid DFT), and C-fitting (for double-hybrid DFT, RI-MP2, DLPNO-CCSD(T)) basis sets. However, as stated above, we are willing to sacrifice some efficiency to have one universal ABS to fit all ERIs. This is something that sets the AutoAux scheme apart from the auto-ABS developed by Yang et al.38 Their ABSs, although large, are designed to fit only Coulomb integrals and are hence applicable only to DFT calculations using pure functionals. Because of the universality requirement, we do not separately examine the errors in the different contributions to the energy but rather use the error in the total

Note that eq 4 is similar to the even-tempered formula (eq 7) described in Section 2. 2. The parameters are then iteratively optimized to minimize either the RI error in the Coulomb (for RIJ)15 or exchange (for RI-JK)17,18 contribution to the energy or the δRI functional (for RI-C),19,22 which is a weighted sum of the squared errors in ERIs that contribute to the MP2 correlation energy. 3. During the optimization process, unnecessary functions are removed, and finally, some functions are contracted. Weigend et al. made a comparison between the CD and RI-X (X = J, JK, or C) approaches in terms of accuracy and efficiency.23 Their conclusion is that the optimized ABSs are sufficiently accurate except for very extended basis sets and do not suffer from the overhead associated with CD procedures. Epifanovsky et al. also benchmarked the performance of the RI and CD approximations applied to coupled-cluster methods.24 They could not identify a clear winner but point out the possibility to fine-tune the accuracy of CD. They also found some issues arising from its dependence on the geometry of the system. In the earlier RI-CCSD(T) implementation, Rendell and Lee simply used the OBS as a fitting basis, noting that this is sufficient but not optimal.25 Recently, Kállay proposed another way to increase the efficiency of RI methods by using tensor decomposition techniques to reduce the size of the ABS, resulting in a system-specific set of so-called natural auxiliary functions (NAFs).26 Both this approach as well as the CD-based ABSs discussed above have the drawback that they are not readily usable in electronic structure codes designed to employ predefined ABSs. The major downside to predefined ABSs is the significant effort required for their design and optimization. The result is that ABSs are often not available or their accuracy has not been ascertained for particular applications or properties. With the development of highly accurate and efficient multireference methodologies,27 the quality of the OBS is increasingly important, with highly augmented or core−valence polarized basis sets being often required. For example, in CASSCF/ NEVPT2 studies of the spectroscopic properties of rare-earth metals, accurate all-electron OBSs are required28 such as the SARC2-QZVP basis sets.29 Although JK-fitting ABSs for the latter have been developed,29 no ABSs exist for similar basis set families such as the Sapporo30−32 or ANO-RCC33−36 basis sets. 555


Article

Journal of Chemical Theory and Computation MP2 energy as a measure of accuracy. While the latter is not variational with respect to the quality of the ABS, its calculation requires all types of ERIs and is therefore sensitive to any deficiencies in the ABS.39 The basic concept of the AutoAux generation procedure derives from the fundamental requirement for an ABS, namely, to span the product space of the OBS as closely as possible. This would be most accurately achieved if there were as many ABFs as there are products of OBFs. Note that the product of two Gaussian basis functions has the form χ1 (r)χ2 (r) =

NY lm1 1(θ ,

ϕ)Y lm2 2(θ ,

−(α1+ α2)r 2

ϕ)exp

l1+ l 2

∑

αl1,

αeff =

2k 2 πrexp

(9)

22l + 1((l + 1) ! )2 (2l + 2)!

(10)

Having obtained a “quasi-uncontracted” OBS, it is a simple matter to examine the orbital product space, find the lowest and highest exponents αlmin and αlmax, respectively, for a given angular momentum l, and construct the even-tempered expansion to encompass them. However, the latter approximation is not entirely justified and results in the auxiliary basis set lacking a few important functions with high exponents, so we multiply the upper limit αlmax by a parameter f l. Finally, to decrease the size of the ABS while keeping the RI error small, we remove some unnecessary ABFs with high angular momenta by applying some heuristic rules, which will be described below and were developed after extensive testing. In the rest of this section, we present a step-by-step algorithm of the AutoAux generation procedure. Given an orbital basis set for element Z: • We find the smallest and largest primitive exponents, l 47 l αmin and αmax,prim , respectively, for every angular momentum l. • We substitute every contracted function with a single Gaussian (with an exponent αeff) which has the same rexp according to eqs 8 and 9. • We then find the largest effective exponent αlmax,eff for every angular momentum l. With this, we describe the boundaries of the OBS and can now work exclusively with the maximal and minimal exponents. • For every pair of angular momenta l and l′ such that l, l′ ∈ [0, lmax] (where lmax is the maximum angular momentum in the OBS) and l ≤ l′: laux = l + l′...|l − aux ′ = αlmin + αlmin , l′|. We store three exponents: αlmin,aux

(5)

(6)

where C(l1, l2, L, m1, m2) are the Clebsch−Gordan (or Wigner) coefficients, which are zero outside the summation range. Therefore, for every pair of OBFs with angular momenta l1 and l2 and exponents α1 and α2, the ABS must include functions with l = |l1 − l2| ... (l1 + l2) and α = α1 + α2. Of course, such an approach would lead to an ABS that is (i) system-dependent, because it includes products of OBFs centered on different atoms, (ii) so large as to defeat the purpose of the RI approximation entirely, and (iii) linearly dependent, because it contains functions with very similar exponents, which leads to numerical instabilities. To resolve the first issue, we disregard products between orbitals on different atomic centers (as is done in the aCD10 and LCAD11 approaches) and can thus produce the ABS starting only from the OBS for a given element. Second, ABFs with the same exponent and angular momentum describe the same electron distribution. Thus, instead of including the actual product functions, we generate an even-tempered basis set spanning the range of exponents given by the product space. An eventempered basis set is defined as having a geometric progression of exponents {αli} for a given angular momentum l: αil = β l αil− 1 , i = 2, ..., N

(8)

k=

C(l1 , l 2 , L , m1 , m2)Y Lm1+ m2

L =|l1− l 2|

∫ χ 2 (r)r dr

where k is a number determined by the angular momentum l of the function and is given by

where r, θ, and ϕ are the polar coordinates, αi is the orbital m exponent, N contains the normalization factors, and Yli i(θ, ϕ) are spherical harmonics dependent on the quantum numbers li and mi. According to the rules for coupling of angular momenta, the product of two spherical harmonics can be expressed as40 Y lm1 1Y lm2 2 =

rexp =

′ aux aux αlmax,aux,prim = αlmax,prim + αlmax,prim , and αlmax,aux,eff = αlmax,eff +

•

(7)

•

with β , and N being the defining parameters. This scheme was first proposed for orbital basis sets by Raffenetti and Ruedenberg,41 and altered versions of it have been employed later,42,43 notably the well-tempered basis sets by Huzinaga and co-workers.44,45 Similar formulations have also been applied to ABSs,14,15,18,46 including eq 4 discussed above. We therefore chose the even-tempered formula for its simplicity and robustness. Another difficulty arises because most Gaussian basis sets are contracted in the core region, and the OBFs with the highest exponents usually have small contraction coefficients. Therefore, spanning the product space of the primitives would result in ABSs that are unnecessarily large. To overcome this issue, we replace every contracted OBF χ with a single Gaussian function (with an exponent αeff) that has the same radial expectation value rexp. To this effect, we employ the formulas: l

•

•

556

′ αlmax,eff . aux and the largest For every laux, we find the smallest αlmin,aux laux laux αmax,aux,prim and αmax,aux,eff. We will use these values to determine the boundaries of the ABS. We assign lval to be the highest angular momentum of the occupied shells in element Z.48 In case an effective core potential (ECP) is employed, we disregard the electrons included in the ECP (e.g., a gold atom with a 60-electron ECP has the electronic configuration 5s25p65d106s1, therefore lval = 2). We limit the maximum angular momentum in the auxiliary basis: lmax,aux = max(2lval, lmax + linc), where linc is 1 for elements up to Ar and 2 for the rest. This is the first heuristic restriction we impose on the ABS. It appears to allow for sufficient accuracy according to our tests but can optionally be lifted by the user. For laux ∈ [0, lmax,aux]: We obtain the actual upper limit aux for the even-tempered expansion αlmax,aux by multiplying laux laux αmax,aux,eff by a factor f > 1 (the actual values used are shown in Table 1). Here, we introduce a second heuristic DOI: 10.1021/acs.jctc.6b01041 J. Chem. Theory Comput. 2017, 13, 554−562

Article


except for MnH, which has a multiplicity of 7 but has been shown to be well-described by single determinant methods.53 The chosen basis sets are def2-SVP, def2-TZVP, def2QZVPP,16,54,55 and cc-pwCVnZ (n = D, T, Q, 5).53,56−66 The respective ECPs are employed for elements beyond Kr (the most recent Stuttgart−Cologne ECPs are used with the ccpwCVnZ basis sets,65−73 while modified versions are used with the def2 basis sets,72−76 as described in ref 55). The geometries of the hydrides were first optimized at the B3LYP-D3/def2TZVP/def2-JK level. Next, single-point MP2 (all-electron and frozen-core) calculations were carried out without RI using the selected OBSs to obtain the total energy, labeled in the following as ENORI. The same calculations were performed with AutoAux using RI-MP2 and RI-JK for the HF part. For comparison, the RI calculations with the def2 basis sets were also done using the def2-JK18 ABS for HF and the appropriate C-fitting ABS (denoted def2-C)22 for MP2. This combination of ABSs is henceforth referred to as def2-JKC.77 The total energies obtained from the RI calculations are denoted EAutoAux and Edef2-JKC, respectively, or ERI in general. The RI error is therefore

Table 1. Default Values of the AutoAux Parameters laux laux

f laux βbig βsmall

0

1

2

3

4

5

6

7

20 1.8

4.0 2.0

4.0 2.2

3.5 2.2

2.5 2.2

2.0 2.3

2.0 3.0

3.0

1.8

approximation. We consider ABFs with laux ≤ 2lval more important in the core region because they potentially describe products of occupied orbitals. On the other hand, ABFs with higher angular momentum allow for polarization in the valence region. We therefore only increase the maximal exponents for shells with low laux while making sure not to overstep the primitive-based aux aux . In short: If laux ≤ 2lval, then αlmax,aux maximum αlmax,aux,prim laux laux laux = max( f αmax,aux,eff, αmax,aux,prim). If laux > 2lval, then aux aux αlmax,aux = αlmax,aux,eff . • Finally, we generate an even-tempered expansion aux and αlNaux ≥ according to eq 7 so that αl1aux = αlmin,aux laux αmax,aux. For the same reason as above, we use a smaller expansion factor for shells with lower angular momentum: if laux ≤ 2lval, we use factor βsmall. Otherwise, we use aux factor βlbig . aux The parameters f laux, βlbig , and βsmall were optimized to reduce the size of the ABS while keeping the RI error in the HF and total MP2 energy (both all-electron and frozen-core) below 200 μEh for the test set of hydrides and noble gases described in the next section.49 The default values in ORCA for these parameters are shown in Table 1.

|ΔE RI| = |E RI − E NORI|

(11)

The absolute value is used because this error can have either sign. The nine reactions in the DC9 benchmark set are described in Table 2 along with their respective reaction energies. Note Table 2. All-Electron MP2 Reaction Energies (in kcal mol−1) for the DC9 Benchmark Seta reaction

3. BENCHMARK CALCULATIONS To assess the accuracy and efficiency of the AutoAux procedure, we performed several sets of benchmark calculations. First, the RI error in the total HF and MP2 (both allelectron and frozen-core) energies was determined using two series of OBSs (the Ahlrichs def2 basis sets and the Peterson cc-pwCVnZ basis sets) of increasing quality for a representative series of elemental hydrides and noble gas atoms, comprising elements from H to Rn. For comparison, the calculations with the def2 OBSs were repeated using the available optimized ABSs. For chemical applications, relative energies are more important than absolute energies. Therefore, we also assessed the accuracy of AutoAux for calculating reaction energies using the DC9 benchmark set.50 Finally, as the main purpose of the RI approximation is to speed up calculations, we measured the computational time required for frozen-core MP2 calculations on cysteine chains of different lengths using OBSs of increasing size and either AutoAux-generated ABSs, optimized ABSs, or no RI approximation. 3.1. Computational Details. All calculations were performed with a development version of the ORCA program package.51 The AutoAux procedure is implemented in the forthcoming ORCA 4.0 release and is available using the keyword “AutoAux” as described in the program manual.52 The set of elemental hydrides and noble gas atoms comprises H2, He, LiH, BeH2, BH3, CH4, NH3, H2O, HF, Ne, NaH, MgH2, AlH3, SiH4, PH3, HS, HCl, Ar, KH, CaH2, MnH, CuH, ZnH2, H2Se, HBr, Kr, RbH, SrH2, AgH, SnH4, HI, Xe, CsH, BaH2, AuH, PbH4, HAt, and Rn. All species are closed-shell

Re1 Re2 Re3 Re4 Re5 Re6 Re7 Re8 Re9

2-pyridone →2-hydroxypyridine (C20)cage → (C20)bowl hepta-1,2,3,5,6-hexaene → hepta -1,3,5-triyne 2 2,3-dimethyl-2-butene → octamethylcyclobutane isomerization of (CH)12 isomerization of carbo-[3]oxacarbon N2CH2 + C2H4 → (CH2)3N2 4Be → Be4 4S2 → S8

def2QZVPP

ccpwCVQZ

ref

−2.8 18.8 −23.4

−2.6 12.6 −23.6

−1.0 −13.3 −14.3

−27.6

−24.3

−19.2

−34.0 −21.7

−30.7 −23.4

−19.5 −26.9

−38.6 −104.6 −184.8

−37.7 −108.4 −182.7

−38.1 −88.4 −101.0

a

Calculated in this work using the def2-QZVPP and cc-pwCVQZ basis sets without the RI approximation and reference values from ref 50.

that the MP2 energies differ significantly from the reference values, most notably for Re2, where even the sign is wrong (as has been described previously78,79). However, in this work, we are not concerned with accurately reproducing these values. Rather, we are interested in the difference between the reaction energies obtained with and without the RI approximation. Therefore, for a model reaction nR → mP (where R and P are the reagents and products and n and m are their respective stoichiometric coefficients), we denote the energies of the reagents and products with ER and EP, respectively, and calculate the value P P R R |ΔΔE RI| = |n(E RI − E NORI ) − m(E RI − E NORI )|

(12)

The energies were once again calculated at the MP2 level of theory (both all-electron and frozen-core) using the def2-SVP, def2-TZVP, def2-QZVPP, and cc-pwCVnZ (n = D, T, or Q) 557


Article

Journal of Chemical Theory and Computation basis sets. For comparison, the reaction energies for the def2 basis sets were also calculated using the def2-JKC basis sets. To assess the efficiency of the generated ABSs, single-point HF and frozen-core MP2 calculations were performed on a cysteine monomer, dimer, and tetramer using the def2-SVP, def2-TZVP, def2-TZVPP, and def2-QZVPP OBSs either with the corresponding optimized ABSs, with AutoAux, or without the RI approximation altogether. The calculations were performed on a single core of an Intel Xeon X5650 2.67 GHz CPU, and the computational times required for the SCF and MP2 parts of the calculation were recorded. For a fairer comparison, enough memory (up to 4 GB for the largest jobs) was provided for the SCF calculations to ensure that all integrals could be calculated in one batch. The geometries of the cysteine chains (in helical conformation) were optimized at the BP86-D3/def2-TZVP/def2-J level. 3.2. Size of the Auxiliary Basis Sets. To estimate the computational efficiency of AutoAux, we generated ABSs to fit each of the examined OBSs as well as a few other commonly used basis sets for all elements where they are available. We then calculated the ratio R = NABS/NOBS for each atom, where N is the number of contracted functions in the respective basis set,80 and took the averages R̅ AutoAux over all atoms with a given OBS. In Table 3, these are compared with the ratios R̅ JK and R̅ C

Figure 1. Number of basis functions per atom of a given element for the def2-TZVP OBS, the corresponding ABS generated by AutoAux, and the optimized JK and C ABSs. Note that def2-TZVP/C basis sets are unavailable for elements Ce−Lu.

AutoAux (see Table 4). Notably, with AutoAux, the average error decreases when increasing the OBS size, whereas it is fairly constant for the optimized ABSs, which is consistent with their design philosophy. It is also worth pointing out that AutoAux performs equally well for all-electron and frozen-core MP2, while the fitted ABSs were meant to be used only with the frozen-core approximation and predictably result in larger errors in the all-electron MP2 energy. A more detailed look at the errors for frozen-core MP2 using the def2 OBSs is provided in Figures 2−4. One can see that the error is usually lower with AutoAux. Notable exceptions include AgH and AuH as well as Kr, SrH2, Pb4, and HAt using def2SVP. However, the error in these cases is still below the threshold value of 200 μEh. The errors for the cc-pwCVnZ basis sets (Figure 5) are also below the threshold with the notable exception of cc-pwCV5Z. As with the def2 basis sets, the RI error in the HF energy decreases for the larger OBSs. However, the error in the MP2 energy increases significantly for cc-pwCV5Z. The most likely reason for this is that for elements beyond Ca the cc-pwCV5Z basis sets include i-functions (l = 6), and the integral package used in ORCA supports functions of angular momentum up to 7, limiting the ABS to lmax + 1 rather than lmax + 2, which would be required to reduce the RI error. In summary, using AutoAux with the parameter set specified in Table 1, almost all of the test cases have RI errors smaller than 200 μEh with average errors in the order of 10−5 Eh, which is as accurate−if not more so−than the optimized ABSs. Errors for the cc-pwCV5Z basis set are higher due to limitations in the highest angular momentum supported by the program. 3.4. Relative Energies: DC9 Benchmark Set. As noted above, because the def2-C basis sets were not meant to be used without the frozen-core approximation, they result in larger RI errors in the all-electron MP2 energies. However, this difference cancels out in the relative energies. Therefore, we only presented results for all-electron MP2 in Table 5 simply because the numbers for frozen-core MP2 are virtually identical. The data show that the errors with AutoAux are similar to the ones with the optimized basis sets: slightly larger for def2-SVP and slightly smaller for def2-TZVP. A breakdown of all the calculated reaction energy errors is given in Figure 6. It is apparent that Re9 (4S2 → S8) is a notable outlier, where AutoAux performs poorly with double-ζ basis sets. A more detailed analysis of the def2-SVP case showed that the error is mainly in the HF energy and cannot be significantly reduced by

Table 3. Average Number of Basis Functions per Atom (N̅ OBS) over All Elements from H to Rn for Which a Given Basis Set is Available and Average Ratios of ABS to OBS Functions for AutoAux and the Fitted JK and C Sets (R̅ AA, R̅ JK, and R̅ C, Respectively) OBS

N̅ OBS

R̅ AA

R̅ JK

R̅ C

def2-SVP def2-TZVP def2-QZVPP cc-pwCVDZc,d cc-pwCVTZc cc-pwCVQZc cc-pwCV5Zc

34 48 87 38 79 129 196

8.41 6.93 5.10 6.15 4.53 3.59 2.74

7.07a 4.80a 2.56a

3.95b 3.36b 2.94b

a

Excluding La, for which the def2-JK basis set is not available. Excluding Ce−Lu, for which def2-C basis sets are not available. c Excluding K, Rb, Cs, and La−Lu, for which cc-pwCVnZ (or ccpwCVnZ-PP, respectively) basis sets are not available. The ccpwCVnZ-PP basis sets are used for elements beyond Kr. dExcluding Sc−Zn, for which cc-pwCVDZ basis sets are not available. b

for the optimized JK and C ABSs, respectively. The results show that AutoAux is up to about twice the size of the optimized ABSs, which is quite large but not prohibitively so and still results in a significant speedup compared to not using RI at all. A more detailed look at the number of basis functions for each element (see Figure 1 for the def2-TZVP case) reveals that this behavior is fairly consistent over the whole periodic table. 3.3. Absolute Energies: Selected Hydrides and Noble Gases. The parameters in the AutoAux procedure were optimized such that the maximal RI errors in the total HF and MP2 (both all-electron and frozen-core) energies for the selected species were less than 200 μEh, which is comparable to using the optimized ABSs. In fact, the average RI error for the HF energy using def2-SVP is about the same with AutoAux as with def2-JKC, while the HF errors with the other OBSs and the MP2 errors with all def2 basis sets are lower when using 558


Article


Table 4. Average RI Error ± Standard Deviation (in μEh) in the Total HF, All-Electron MP2 (AE-MP2), and Frozen-Core MP2 (FC-MP2) Energies Using AutoAux or the Optimized ABSs (def2-JKC) AutoAux

def2-JKC

OBS

HF

FC-MP2

AE-MP2

HF

FC-MP2

AE-MP2

def2-SVP def2-TZVP def2-QZVPP cc-pwCVDZ cc-pwCVTZ cc-pwCVQZ cc-pwCV5Z

28 ± 37 8±9 3 ± 10 32 ± 52 4±7 2±5 2±5

53 ± 45 25 ± 24 17 ± 22 53 ± 55 29 ± 38 44 ± 50 125 ± 178

50 ± 42 24 ± 23 17 ± 21 45 ± 38 32 ± 43 50 ± 51 175 ± 194

33 ± 25 22 ± 20 24 ± 23

89 ± 56 86 ± 51 87 ± 69

100 ± 56 95 ± 50 162 ± 138

Figure 2. Unsigned RI error in the total frozen-core MP2 energy of the selected hydrides, calculated using the def2-SVP basis set.

Figure 3. Unsigned RI error in the total frozen-core MP2 energy of the selected hydrides, calculated using the def2-TZVP basis set.

error of 0.8 kcal mol−1 even less significant. Because of these considerations, we deem the larger error acceptable. The same reasoning applies to the cc-pwCVDZ calculation as well. The situation is somewhat less clear for the Re4 calculation with def2-QZVPP/AutoAux, where removing some functions from the ABS actually decreases the RI error, suggesting that it is due to unfortunate and unsystematic addition of errors. In short, the RI errors in relative energies are under 1 kcal mol−1 (by an order of magnitude in most cases), which is sufficiently small for chemical applications and negligible in comparison with the other errors (due to the method, basis set incompleteness and superposition, etc.) associated with the calculations.

changing the AutoAux parameters. It was found that adding the most diffuse functions from the def2-JK ABS to the generated one greatly reduces the RI error. However, these functions are outside the range of the atomic basis product space, so there is no general way to include them in the generated ABS without taking the molecular geometry into account. Therefore, we suggest caution when using AutoAux with small basis sets. Note, however, that the lack of diffuse functions in the def2SVP basis set also results in a large basis set error: the reaction energy is calculated to be −140.4, −178.0, and −184.8 kcal mol−1 with the def2-SVP, def2-TZVPP, and def2-QZVPP OBSs, respectively (the reference value is −101.0 kcal mol−1). Furthermore, this large reaction energy makes the already small 559


Article


Figure 4. Unsigned RI error in the total frozen-core MP2 energy of the selected hydrides, calculated using the def2-QZVPP basis set.

Figure 5. Unsigned RI error in the total all-electron MP2 energy of the selected hydrides, calculated using the cc-pwCVnZ (n = D, T, Q, or 5) basis sets.

Table 5. Average RI Error ± Standard Deviation (in kcal mol−1) in the All-Electron MP2 Reaction Energies for the DC9 Benchmark Set Using AutoAux or the Optimized ABSs (def2-JKC) OBS def2-SVP def2-TZVP def2-QZVPP cc-pwCVDZ cc-pwCVTZ cc-pwCVQZ

AutoAux

def2-JKC

± ± ± ± ± ±

0.06 ± 0.03 0.05 ± 0.05 0.05 ± 0.06

0.16 0.03 0.05 0.10 0.03 0.01

0.23 0.03 0.11 0.13 0.03 0.01

3.5. Performance: Timings for Cysteine Chains. While the calculations discussed in Section 3.3 are fast enough to be feasible without the RI approximation, it is important to know what speedup can be achieved using RI in conjunction with AutoAux in a more realistic scenario and how this compares to using optimized ABSs. For this purpose, timings for HF and MP2 calculations on cysteine chains were obtained and are presented in Table 6. All SCF calculations converged in 10−14 iterations with the RI calculations typically requiring 1−3 cycles more than the analytic ones. The RI-MP2 calculations with AutoAux are an order of magnitude faster than the conventional ones and up to about 2.5 times slower than the ones with

Figure 6. Unsigned RI errors in the all-electron MP2 reaction energies for the DC9 benchmark set. 560


Article


Table 6. Total CPU Time (in Seconds) Required to Calculate the HF and MP2 Parts of Frozen-Core MP2 Calculations on Cysteine Chains of Different Lengths Using Either No RI or RI-JK and RI-MP2 with AutoAux or the Corresponding Optimized ABSs def2-SVP (Cys)1

(Cys)2

(Cys)4

a

no RI AutoAux def2-JKC no RI AutoAux def2-JKC no RI AutoAux def2-JKC

def2-TZVP

def2-TZVPP

def2-QZVPP

HF

MP2

HF

MP2

HF

MP2

HF

MP2

38 36 27 178 261 175 1076 2898 1841

15 2 1 101 22 9 1207 362 154

371 162 66 2046 1313 474 11 852 17 377 5042

167 12 5 1262 143 55 14 758 2566 1047

607 237 86 3168 1944 619 17 602 24 506 6145

297 19 8 2231 233 95 27 320 4264 1815

6834 1297 321 38 456 9236 1951 a

3735 101 52 24 783 1196 623

Calculations were not performed on (Cys)4 with the def2-QZVPP basis set.

MP2 and other correlated methods. Our goal was to provide users of ORCA with a robust and universal tool that can be used when no auxiliary basis set is available. In this we believe to have been successful: we showed that in reproducing absolute and relative MP2 energies, AutoAux is as accurate or more so (with the few exceptions discussed above) than the optimized ABSs for the def2 series of basis sets and also works well for the cc-pwCVnZ basis sets, which include tight core polarization functions. The simplicity of the scheme comes at the cost that the generated ABSs can be about twice the size of the optimized ones but are still small enough to benefit from use of the RI approximation, especially for MP2 calculations. Nevertheless, it is possible to improve on the efficiency of the procedure, and we intend to refine it in future work.

the def2-C ABSs, which is to be expected considering the larger size of AutoAux (the AutoAux ABSs for the def2-SVP, def2TZVP, def2-TZVPP, and def2-QZVPP OBSs are respectively 1.35, 2.19, 2.44, and 3.52 times larger than the def2-JK ABS and 2.15, 2.31, 2.21, and 1.79 times larger than the corresponding def2-C ABSs). The situation is somewhat more complicated for the HF part: it has already been established that the improvement of RI-JK over conventional HF is highest for small to midsize systems with large (over triple-ζ) basis sets.17,18 In the given case, the def2-QZVPP/def2-JK calculations are up to about 20 times faster than the conventional ones, while more modest factors of 2−5 are obtained for the larger systems and smaller OBSs. Note that the def2-SVP calculation on (Cys)4 takes longer with RI than without. This is to be expected, as the formal scaling of RI-JK is still N4, and the efficiency increases proportionally to the ratios NOBS/Nocc (Nocc being the number of occupied orbitals) and NOBS/NABS.17 As a consequence, this crossover occurs earlier with AutoAux due to the larger size of the ABS, so a longer computational time compared to that of the standard procedure is required for all calculations on (Cys)4 as well as for (Cys)2 with the def2-SVP basis set. While this does limit the usefulness of AutoAux for RI-JK HF calculations to smaller systems than would be feasible with the optimized ABSs, this is not a limitation of AutoAux as such but is simply a result of the scaling of the RI-JK procedure. For larger systems, one would be well-advised to instead employ methods with superior scaling, such as RIJCOSX.81 For example, the RIJCOSX calculations on (Cys)4 with the def2-SVP, def2-TZVP, and def2-QZVPP basis sets took, respectively, 944, 4242, and 5137 s with the def2-JK basis set and 951, 4411, and 5383 s with AutoAux,82 showing only a minor dependence on the ABS. In summary, RI-MP2 calculations using AutoAux are up to 2.5 times slower than those using the optimized basis sets but are still significantly faster than exact MP2 calculations. For HF, the use of RI-JK with AutoAux is beneficial only for fairly small systems (less than 200 electrons) with large basis sets (at least triple-ζ quality). For bigger systems, RIJCOSX is recommended over RI-JK.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Georgi L. Stoychev: 0000-0002-3408-5749 Funding

The authors are grateful to the Max Planck Society for generous support of our work. G.L.S. thanks the European Commission for financial support through the Erasmus+ program. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Almlöf, J.; Faegri, K.; Korsell, K. J. Comput. Chem. 1982, 3, 385− 399. (2) Häser, M.; Ahlrichs, R. J. Comput. Chem. 1989, 10, 104−111. (3) Almlöf, J.; Faegri, K. In Self-Consistent Field: Theory and Applications; Carbó, R., Klobukowski, M., Eds.; Elsevier: Amsterdam, 1990; p 195. (4) Horn, H.; Weiß, H.; Háser, M.; Ehrig, M.; Ahlrichs, R. J. Comput. Chem. 1991, 12, 1058−1064. (5) Almlöf, J. In Modern Electronic Structure Theory; Yarkony, D. R., Ed.; World Scientific: Singapore, 1995; Chapter 3, pp 110−151. (6) Whitten, J. L. J. Chem. Phys. 1973, 58, 4496. (7) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, 71, 3396. (8) Vahtras, O.; Almlöf, J.; Feyereisen, M. W. Chem. Phys. Lett. 1993, 213, 514−518. (9) Beebe, N. H. F.; Linderberg, J. Int. J. Quantum Chem. 1977, 12, 683−705.

4. CONCLUSION The presented work describes the AutoAux procedure, which has been implemented in the forthcoming version 4.0 of the ORCA computational chemistry package and generates an auxiliary basis set on-the-fly to be used with the RI approximation of ERIs appearing in HF and DFT as well as 561


Article


polarization functions) to make sure we fit the valence shell as well as possible. (48) After testing, we found it necessary to assign lval = 1 for Li and Be. (49) The value of 200 μEh was chosen so that the errors with AutoAux and the optimized basis sets would be comparable. Note, however, that because the total energy is an extensive quantity, the RI error is expected to increase with the system size. (50) Goerigk, L.; Grimme, S. J. Chem. Theory Comput. 2010, 6, 107− 126. (51) Neese, F. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 73− 78. (52) ORCA User Manual. https://orcaforum.cec.mpg.de/ OrcaManual.pdf. To be released with ORCA version 4.0. (53) Balabanov, N. B.; Peterson, K. A. J. Chem. Phys. 2005, 123, 064107. (54) Weigend, F.; Furche, F.; Ahlrichs, R. J. Chem. Phys. 2003, 119, 12753−12762. (55) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297. (56) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (57) Woon, D. E.; Dunning, T. H. J. Chem. Phys. 1994, 100, 2975. (58) Koput, J.; Peterson, K. A. J. Phys. Chem. A 2002, 106, 9595− 9599. (59) Peterson, K. A.; Dunning, T. H. J. Chem. Phys. 2002, 117, 10548−10560. (60) Peterson, K. A.; Puzzarini, C. Theor. Chem. Acc. 2005, 114, 283− 296. (61) DeYonker, N. J.; Peterson, K. A.; Wilson, A. K. J. Phys. Chem. A 2007, 111, 11383−11393. (62) Peterson, K. A.; Yousaf, K. E. J. Chem. Phys. 2010, 133, 174116. (63) Prascher, B. P.; Woon, D. E.; Peterson, K. A.; Dunning, T. H.; Wilson, A. K. Theor. Chem. Acc. 2011, 128, 69−82. (64) Li, H.; Feng, H.; Sun, W.; Zhang, Y.; Fan, Q.; Peterson, K. A.; Xie, Y.; Schaefer, H. F. Mol. Phys. 2013, 111, 2292−2298. (65) Peterson, K. A.; Figgen, D.; Dolg, M.; Stoll, H. J. Chem. Phys. 2007, 126, 124101. (66) Figgen, D.; Peterson, K. A.; Dolg, M.; Stoll, H. J. Chem. Phys. 2009, 130, 164108. (67) Metz, B.; Schweizer, M.; Stoll, H.; Dolg, M.; Liu, W. Theor. Chem. Acc. 2000, 104, 22−28. (68) Lim, I. S.; Schwerdtfeger, P.; Metz, B.; Stoll, H. J. Chem. Phys. 2005, 122, 104103. (69) Figgen, D.; Rauhut, G.; Dolg, M.; Stoll, H. Chem. Phys. 2005, 311, 227−244. (70) Peterson, K. A.; Shepler, B. C.; Figgen, D.; Stoll, H. J. Phys. Chem. A 2006, 110, 13877−13883. (71) Lim, I. S.; Stoll, H.; Schwerdtfeger, P. J. Chem. Phys. 2006, 124, 034107. (72) Metz, B.; Stoll, H.; Dolg, M. J. Chem. Phys. 2000, 113, 2563. (73) Peterson, K. A.; Figgen, D.; Goll, E.; Stoll, H.; Dolg, M. J. Chem. Phys. 2003, 119, 11113. (74) Andrae, D.; Häußermann, U.; Dolg, M.; Stoll, H.; Preuß, H. Theor. Chem. Acc. 1990, 77, 123−141. (75) Kaupp, M.; Schleyer, P. v. R.; Stoll, H.; Preuss, H. J. Chem. Phys. 1991, 94, 1360. (76) Leininger, T.; Nicklass, A.; Küchle, W.; Stoll, H.; Dolg, M.; Bergner, A. Chem. Phys. Lett. 1996, 255, 274−280. (77) Note that the def2-C basis sets are different for def2-SVP, def2TZVP, and def2-QZVPP, unlike the def2-JK sets, which are universal. (78) Taylor, P. R.; Bylaska, E.; Weare, J. H.; Kawai, R. Chem. Phys. Lett. 1995, 235, 558−563. (79) An, W.; Gao, Y.; Bulusu, S.; Zeng, X. C. J. Chem. Phys. 2005, 122, 204109. (80) Note that we are counting functions with different values of m separately. (81) Kossmann, S.; Neese, F. Chem. Phys. Lett. 2009, 481, 240−243. (82) A 194-point Lebedev final grid was used to ensure comparable errors in energy; see the discussion in ref 13 for details on grid sizes.

(10) Aquilante, F.; Lindh, R.; Bondo Pedersen, T. J. Chem. Phys. 2007, 127, 114107. (11) Ten-no, S.; Iwata, S. Chem. Phys. Lett. 1995, 240, 578−584. (12) Ten-no, S.; Iwata, S. Int. J. Quantum Chem. 1996, 60, 1319. (13) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Chem. Phys. 2009, 356, 98−109. (14) Eichkorn, K.; Treutler, O.; Ö hm, H.; Häser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 240, 283−290. (15) Weigend, F. Phys. Chem. Chem. Phys. 2006, 8, 1057−65. (16) Gulde, R.; Pollak, P.; Weigend, F. J. Chem. Theory Comput. 2012, 8, 4062−4068. (17) Weigend, F. Phys. Chem. Chem. Phys. 2002, 4, 4285−4291. (18) Weigend, F. J. Comput. Chem. 2008, 29, 167−175. (19) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294, 143−152. (20) Weigend, F.; Köhn, A.; Hättig, C. J. Chem. Phys. 2002, 116, 3175−3183. (21) Hättig, C. Phys. Chem. Chem. Phys. 2005, 7, 59−66. (22) Hellweg, A.; Hättig, C.; Höfener, S.; Klopper, W. Theor. Chem. Acc. 2007, 117, 587−597. (23) Weigend, F.; Kattannek, M.; Ahlrichs, R. J. Chem. Phys. 2009, 130, 164106. (24) Epifanovsky, E.; Zuev, D.; Feng, X.; Khistyaev, K.; Shao, Y.; Krylov, A. I. J. Chem. Phys. 2013, 139, 134105. (25) Rendell, A. P.; Lee, T. J. J. Chem. Phys. 1994, 101, 400−408. (26) Kállay, M. J. Chem. Phys. 2014, 141, 244113. (27) Guo, Y.; Sivalingam, K.; Valeev, E. F.; Neese, F. J. Chem. Phys. 2016, 144, 094111. (28) Aravena, D.; Atanasov, M.; Neese, F. Inorg. Chem. 2016, 55, 4457−4469. (29) Aravena, D.; Neese, F.; Pantazis, D. A. J. Chem. Theory Comput. 2016, 12, 1148−1156. (30) Noro, T.; Sekiya, M.; Koga, T. Theor. Chem. Acc. 2012, 131, 1124. (31) Sekiya, M.; Noro, T.; Koga, T.; Shimazaki, T. Theor. Chem. Acc. 2012, 131, 1247. (32) Noro, T.; Sekiya, M.; Koga, T. Theor. Chem. Acc. 2013, 132, 1363. (33) Widmark, P.-O.; Malmqvist, P.-Å; Roos, B. O. Theor. Chem. Acc. 1990, 77, 291−306. (34) Roos, B. O.; Lindh, R.; Malmqvist, P.-Å; Veryazov, V.; Widmark, P.-O. J. Phys. Chem. A 2005, 109, 6575−6579. (35) Roos, B. O.; Lindh, R.; Malmqvist, P.-Å.; Veryazov, V.; Widmark, P.-O. J. Phys. Chem. A 2004, 108, 2851−2858. (36) Roos, B. O.; Lindh, R.; Malmqvist, P.-Å.; Veryazov, V.; Widmark, P.-O. Chem. Phys. Lett. 2005, 409, 295−299. (37) Rappoport, D.; Furche, F. J. Chem. Phys. 2010, 133, 134105. (38) Yang, R.; Rendell, A. P.; Frisch, M. J. J. Chem. Phys. 2007, 127, 074102. (39) During our work, we also separately monitored the HF energy and MP2 correlation energy to not be mislead by fortuitous error cancellation. Only fitting the Coulomb part (RI-J) using AutoAux also leads to similar errors in the HF energy. These results have been omitted here for the sake of brevity. (40) Rose, M. E. Elementary Theory of Angular Momentum; John Wiley & Sons: New York, 1957; p 272. (41) Raffenetti, R.; Ruedenberg, K. Even-tempered representations of atomic self-consistent-field wavefunctions, Report IS-3195; USAEC, Iowa State University: Ames, IA, 1973. (42) Bardo, R. D.; Ruedenberg, K. J. Chem. Phys. 1974, 60, 918. (43) Clementi, E.; Corongiu, G. Chem. Phys. Lett. 1982, 90, 359−363. (44) Huzinaga, S.; Miguel, B. Chem. Phys. Lett. 1990, 175, 289−291. (45) Huzinaga, S.; Klobukowski, M. Chem. Phys. Lett. 1993, 212, 260−264. (46) Früchtl, H. A.; Kendall, R. A.; Harrison, R. J.; Dyall, K. G. Int. J. Quantum Chem. 1997, 64, 63−69. (47) Note the lack of an index “prim”. This is because we always start the even-tempered expansion from the smallest primitive exponent, even if it is part of a contraction (which is rarely the case except for 562


Automatic Generation of Auxiliary Basis Sets - ACS Publications

Recommend Documents