Gaussian Basis Sets for Crystalline Solids: All-Purpose Basis Set Libraries vs System-Specific Optimizations

It is customary in molecular quantum chemistry to adopt basis set libraries in which the basis sets are classified according to either their size (triple-ζ, quadruple-ζ, ...) and the method/property they are optimal for (correlation-consistent, linear-response, ...) but not according to the chemistry of the system to be studied. In fact the vast majority of molecules is quite homogeneous in terms of density (i.e., atomic distances) and types of bond involved (covalent or dispersive). The situation is not the same for solids, in which the same chemical element can be found having metallic, ionic, covalent, or dispersively bound character in different crystalline forms or compounds, with different packings. This situation calls for a different approach to the choice of basis sets, namely a system-specific optimization of the basis set that requires a practical algorithm that could be used on a routine basis. In this work we develop a basis set optimization method based on an algorithm–similar to the direct inversion in the iterative subspace–that we name BDIIS. The total energy of the system is minimized together with the condition number of the overlap matrix as proposed by VandeVondele et al. [VandeVondele et al. J. Chem. Phys.2007, 227, 114105]. The details of the method are here presented, and its performance in optimizing valence orbitals is shown. As demonstrative systems we consider simple prototypical solids such as diamond, graphene sodium chloride, and LiH, and we show how basis set optimizations have certain advantages also toward the use of large (quadruple-ζ) basis sets in solids, both at the DFT and Hartree–Fock level.


INTRODUCTION
When dealing with the quantum chemical modeling of crystalline solids, the existence of various types of chemical bonding is clearly evident. For instance, the polymorphism of carbon in the graphite (or graphene) and diamond allotropes is just one of many examples, in which the profoundly different chemical behavior is manifested by the same chemical element in different crystal packings. Another exemplary case is that of rocksalt NaCl: sodium is by nature metallic as a bulk material, and chlorine is commonly found in the form of a molecular crystal Cl 2 . NaCl is a prototypical ionic salt. The chemical differences in those materials can be made evident by looking at their electron density (see Figure 1): the electrons involved in the metallic bond in Na are quite spread out over the whole space, while in Cl 2 the density is somewhat more localized on molecules, with empty space between them. Conversely, the wave function in an ionic system like NaCl is strongly confined in a vicinity of the ions and features nodes in the planes in between neighboring atoms. NaCl is also considerably more densely packed.
This variety of chemical bondings in the solid state then reflects the choice of the type and quality of the basis set adopted in the mathematical form of the wave function when solving the Schrodinger equation within periodic boundary conditions (i.e., Bloch functions). 1−3 The situation in the field of molecular modeling is somewhat simpler as isolated molecules or molecular aggregates have nearly comparable atomic densities, and there are commonly no analogue extended systems featuring metallic, ionic, or covalent bonds. Therefore, in molecular calculations, atom-centered basis sets as Gaussian-type orbitals 4 are almost universally adopted, 5 although other basis sets can be and are eventually used.
On the other hand, for solid-state calculations, 2 plane waves, 6−8 atom-centered Gaussians 9 (or their combinations 10 ), and numerical basis sets 11,12 are all popular choices. The plane wave basis, that is naturally suited for nonlocal wave functions such as in the uniform electron gas or in a metal, has the undeniable advantage of a one-knob tuning of accuracy and cost through the kinetic energy cutoff parameter. However, the correct description of local orbitals, core states, or the void can result in a rather high computational cost. Similarly, the inclusion of exact HF exchange in hybrid HF/DFT calculations leads to a steep increase in computational time. Gaussian-type basis sets are less commonly adopted for the quantum chemical treatment of solids, with respect to plane waves. Gaussian functions have the great advantage of allowing to transfer to the solid state a large part of the technology and knowledge that is the legacy of several decades of advances in molecular quantum chemistry and to retain the chemical intuition when looking at the electronic charge distribution of the investigated system. The price to pay is the mandatory definition of a basis set for each atomic species, that is ultimately left in the hands of the end user.
Nowadays, standardized basis set libraries are not commonly available for solids as they are for molecules, 13,14 despite recent attempts in that direction being carried out by Bredow and coworkers. 15−17 The reasons are not only to be ascribed to a lesser effort in a systematic construction of all-purpose basis sets but also more specifically to the wide difference in chemical bonding as outlined above. First attempts to understand the role of basis functions in solids were done by Hess and co-workers, 18 but also more recently Jensen 19 compared atomic, molecular, and solid-state basis sets for carbon and silicon to highlight the differences originating from the different chemical environments.
Another aspect related to the adoption of Gaussian-type functions is the basis set incompleteness due to the use of a finite number of basis functions. Basis set incompleteness is an issue in all types of calculations, but most of all in calculations that employ atom-centered basis sets−Gaussians, Slater functions, or numerical orbitals. This is because the atomic basis sets can never be made complete enough in polyatomic systems, as the basis becomes overcomplete−necessitating the removal of variational degrees of freedom−before becoming complete. In molecules it is rather common to adopt a sequence of basis sets of increasing size (e.g., cc-pVXZ (X = D,T,Q,···) 20 and pc-X (X = 1,2,3, ...) 21 ), but this is not yet routinely applicable for solids. Therefore, reaching the basis set limit is not trivial−even for such simple systems as lithium hydride 22−26 −and is not just a matter of computational efforts: as basis sets grow larger, exponents tends to become more diffuse, linear dependency problems arise, and the convergence of infinite Coulomb and exchange series is jeopardized.
The problem of linear dependencies with an extended basis set is a matter of active research not only for solids but also for average-sized molecules. 27 While the important role of diffuse functions in solids has been recently highlighted by Kadek et al., 28 too diffuse functions are often not needed for ground state calculations because of the packing of the atoms in the unit cell. Such very diffuse functions can also be added a posteriori through dual basis set techniques. 29 Seen from another viewpoint, the main conceptual difference in basis sets meant for the solid state as opposite to molecular electronic structure calculations is that the latter have to describe the asymptotic exponential decay of the electron density in a finite system, requiring somewhat diffuse functions, whereas diffuse basis functions are generally thought not to be necessary in solid-state calculations because the density is much more uniform throughout the cell. In this work our aim is to (i) show to what extent the basis sets are different in different chemical environments, by optimizing bases of the def2-TZVP quality 30−32 and (ii) attempt to use suitably optimized quadruple-ζ basis sets, also from the def2-family, to verify whether they can be adopted for solids without significant pruning, and outline possible strategies for reaching such goal.
To this purpose we present a technique for the optimization of basis set exponents and contraction coefficients, that is based on the Direct Inversion in the Iterative Subspace (DIIS) technique 33−35 and actually quite similar to its geometry optimization variant, GDIIS. 36 The algorithm is implemented in the CRYSTAL code. 9 We show how such optimization allows the retaining of the full number of Gaussians letting the algorithm decide about the diffuseness of the exponents.

THEORETICAL FRAMEWORK
2.1. Background. In the linear combinations of atomic orbitals (LCAO) framework, the crystalline orbitals ψ are treated as linear combinations of Bloch functions (BF) ϕ that are, in turn, defined in terms of local atom-centered functions in which g is a direct space lattice vector, k is the lattice vector defining a point in the reciprocal lattice, A are the coordinates of the atom in the reference cell on which the AO φ is centered, and a are the variational coefficients. The sum over μ is limited to the number of basis functions in the unit cell. The Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article sum over g is, in principle, extended to all the (infinite) lattice vectors of the direct lattice; therefore, suitable screening techniques have to be adopted. 1,38,39 As usual, the AOs can be written as a contraction of a number of primitive Gaussian-Type Functions (GTF) G centered on the same atom in which d j are the contraction coefficients, and α j are the exponents of the radial component of the function. The number, type, and contraction scheme of the Gaussian basis set define its quality. Gaussian functions are defined as where R l (r) = Ne −α·r 2 is the radial part−N being a normalization constant−and Y lm (θ, ρ) is a spherical harmonic.
2.2. The BDIIS Method. Our goal is to devise a suitable algorithm for a system-specific optimization of the exponents α j and contraction coefficients d j as in eq 3. Taking inspiration from the well-known Direct Inversion of Iterative Subspace (DIIS) algorithm of Pulay, 33,34 we describe in the following our Basis-set DIIS (BDIIS) method.
The idea is that of an iterative procedure in which, at each step n, exponents and contraction coefficients are obtained as a linear combination of the trial vectors obtained in previous iterations In the above, e i α and e i d are, respectively, the changes in exponents and contraction coefficients, as predicted by a simple Newton−Raphson step. In fact the gradients e i are defined by where Ω is a suitable functional to be minimized. Here we decide to minimize the system's total energy to which we add a penalty function including the Overlap matrix condition number, following the proposal of VandeVondele and Hutter: 40 The value γ = 0.001 was adopted as suggested in ref 40. In eq 8, κ({α, d}) is the condition number, i.e., the ratio between the largest and the smallest eigenvalue of the overlap matrix at the center of the Brillouin zone (Γ-point). The purpose of such penalty function is to prevent the onset of harmful linear dependence. Linear dependence issues can give rise to numerical instabilities and, as a consequence of that, the appearance of unphysical states. Such unphysical states generally lead to a catastrophic behavior of the total energy that can drop to a value that is orders of magnitude larger, in absolute value, than the proper one.
Although the first of derivatives in (7) could be in principle computed analytically, 41,42 in the present work we evaluate both e i α and e i d by means of numerical derivatives (vide infra). The length of the estimated Newton steps represented by the e α and e d can assume the meaning of an estimated distance from the minimum of Ω and thus be utilized as a measure of the "error" at step n.
The DIIS error matrix, that has the size of the iterative space considered, is built from the scalar products , we can obtain the linear combination coefficients of the BDIIS method to be used in (5) and (6) by solving the linear equation system i k j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j y where λ is a Lagrange multiplier. Such an approach is, in fact, similar to geometry-optimization DIIS (GDIIS 36 ) adopting an identity Hessian. 2.3. Details of the Implementation. The BDIIS procedure outline above has been implemented in a development version of the CRYSTAL17 code. 9 As already mentioned, in this work we compute the derivatives in (7) by means of a twosided numerical derivative. Which means, for exponents α (11) and similarly for coefficients d.
We have also tried to compute a diagonal Hessian using the three points α i + Δα ̅ , α i and α i − Δα ̅ , so to improve the step (error) as defined in eqs 5 and 6 at the same computational cost. However, such a diagonal Hessian seemed not to improve on the quality of the step, and the overall convergence pattern turned out to be similar or slower in all cases we tested. We surmise that the cause can reside in the insufficient accuracy of a three-point numerical estimate of the second derivative.
Once a suitable step Δα ̅ n = α ̅ n − α n−1 is obtained from eq 5, a line search is performed for tuning the optimal parameter f l f n n l n 1 by sampling f l from 0.1 to 1 in a suitable discrete point grid.
The point with the minimum value of Ω is then retained. The convergence of the iterative optimization procedure is verified by checking the absolute value of the largest component of both the gradients and the penalty function. The iterative space used in the BDIIS procedure is set at most to the 14 previous cycles, and the BDIIS step is active since the second basis set optimization step. The optimization is complete when the absolute value of the difference in the penalty function is less than 1.0 · 10 −5 au and the absolute value of the largest component of gradient converges to 3.0 · 10 −4 .

RESULTS
In this section we first briefly describe the performance of the BDIIS method in minimizing the Ω energy functional as Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article defined in eq 8. Then we focus on the effect of system-specific basis set optimizations, by showing the differences between optimized exponents of a typical triple-ζ basis in simple systems containing the same atoms but in a different chemical bonding situation. Finally, we analyze how extended basis sets, such as molecular quadruple-ζ quality, can be optimized for dense solids without significant pruning. In the Supporting Information the reader can find full CRYSTAL17 9 inputs for all the calculations presented in the following, including the explicit definition of the atomic basis sets. All of the optimizations in this work have been carried out starting from molecular def2-TZVP or def2-QZVP basis sets. 30−32 Although the implemented algorithm is general, as described in the previous section, in the following we will focus on the optimization of valence and polarization functions only−the ones relevantly changing in a different chemical environment. Since they are usually uncontracted Gaussian functions, the optimization has been performed solely for the exponents. As a general strategy, we took as a starting point the molecular basis sets, upscaled the exponents of all outermost functions so to avoid small values (<0.1) but without pruning the basis set, and finally optimized the corresponding values by minimizing the function Ω.
In many cases, we adopted pure GGA functionals such as PBE 43 and PBEsol 44 in order to have a faster time to solution. In other cases we used PBE0 45 or Hartree−Fock. More generally, we do not regard our basis set optimizations to be much dependent on the chosen method, 46 since we do not deal with the reoptimization of the core. As the focus of our work is on accuracy and numerical stability, we will not present timings.
3.1. Performance of the BDIIS Method. In Figure 2 we report the progress of the Ω functional minimization−cf. eq 8−along with the BDIIS iterations, in two exemplary yet challenging cases for Gaussian-type basis sets: graphene and bulk metallic sodium. In graphene, the basis set optimizer, run with the PBEsol functional, leads to a stable result after a few iterations, which represents a significant energy gain with respect to the starting point and remains stable for long. If the optimization is allowed to continue for hundreds of cycles, a rise in the penalty function γ ln κ({α,d}) is observed, which evidently prevents the gaussians to become too diffuse. A corresponding decrease of the electronic energy is observed. We remark that such changes are however minimal with respect to the effect of the first iterations, and the optimization is essentially converged after 50 cycles to all practical purposes. In the same figure we have also reported the curve obtained using the Broyden−Fletcher−Goldfarb−Shanno (BFGS) method. It is seen that such a method reaches the same value of the Ω functional, more slowly but also more stably. We will discuss the differences in the solution in the following.
The case of bulk metallic sodium (Figure 2(a)) is different: the electronic energy varies little (and even increases slightly with respect to the starting point); but the penalty function is much more relevant than in other cases, and about 100 iterations are required to reach a plateau. Notably, in this case the basis set optimization was carried out with a hybrid HF/ DFT functional (i.e., PBE0). This level of theory is usually expected to be problematic for metallic systems, but the BDIIS algorithm runs smoothly to convergence.
3.2. Role of the Chemical Environment. We compare here two sets of systems, composed by the same elements: first crystalline diamond, graphene, and carbyne chain and then NaCl with bulk Na and Cl solids. We compare our systemdependent optimized basis sets with the pob-TZVP 15,17 ones. These were also derived from def2-TZVP but differently from ours: (i) the valence exponents were optimized for each system in a comprehensive set of solids with different chemical environments, (ii) for multiple optimization of the same atomic species an averaged value of the exponent was considered, and (iii) most notably, many of the outermost functions were removed, thus reducing the consistent quality of the basis.
We will refer to the basis sets optimized in this work as "dcm-TZVP". Since different basis sets of the same nominal Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article quality are obtained through optimization on different systems, we will adopt the more detailed notation dcm[···]-TZVP, specifying in square brackets the system used for the optimization (e.g., dcm[NaCl]-TZVP).
3.2.1. Diamond, Graphene, and Carbyne Chain. Diamond and graphene are two allotropes of carbon. Both are covalently bound systems but differ by hybridization (sp 3 and sp 2 ), as well as crystalline (3D Vs 2D) and electronic (insulator and conductor) structures. Carbyne is a model system with 1D periodicity (polymer), two atoms in the unit cell and alternating bond length.
In Table 1 we compare the exponents of the original def2-TZVP, the original pob-TZVP, the recently revised pob-rev2 basis set, and our dcm-TZVP basis specifically optimized for diamond, graphene, and carbyne with the PBE functional. For brevity, we will refer to the latter two basis sets as dcm[C diam ]-TZVP, dcm[C graph ]-TZVP, and dcm[C cby ]-TZVP, respectively. Figure 3 shows a corresponding graphical representation of the radial component of some of the involved gaussians. The first striking effect observed is the overall contraction of exponents with respect to the molecular basis. This is not unexpected 15,19 and is to be ascribed to the higher density of atoms in the solid-state phase.
The outermost p-type function shows probably the most significant difference between diamond and graphene. Such difference is due both to the different chemical bonding (sphybridization) and atomic density−graphene is a 2D system surrounded by vacuum in the third dimension. This vacuum offers more space for the Gaussian functions to expand and at the same time requires more extended functions to cover that empty space. Such interpretation is corroborated by the example of the 1D carbyne chain basis dcm[C cby ]-TZVP, which features an even more diffuse p-shell. We take the opportunity here to remind that−conversely to plane-waves−in an atomcentered Gaussian-based approach the true 2D and 1D periodicity is possible, hence the vacuum in the nonperiodic directions is a true vacuum. The effect of the reduced dimensionality is, however, partly counterbalanced by a progressively shorter carbon−carbon distance that is 2.92 Å in diamond, 2.69 Å in graphene, and 2.39/2.46 Å in carbyne due to the different hybridization of the carbon atom in the three compounds. The more diffuse p-function is responsible for the failed convergence when using the graphene dcm-[C graph ]-TZVP basis set in diamond (Table 2). Also d-and ftype functions have a somewhat different spread in the two systems, showing that quadrupole and octupole interactions act differently in the two allotropes.
In Table 2 we report some total energies obtained at the DFT/PBE level: in addition to dcm-TZVP and pob-TZVP bases, the dcm[C diam ]-TZVP basis was also tested in graphene and the dcm[C graph ]-TZVP in diamond. From Table 2, we see that the energies relative to the proper dcm bases are lower by about 0.014 E h than the pob-ones. On the other hand, swapping the two dcm-TZVP bases led to an energy similar to (though still lower than) that of pob-TZVP[G] while the more diffuse dcm[C graph ]-TZVP turned out to be unusable in the   Let us now compare the optimal basis set obtained with the PBE0 functional for three bulk structures with very different chemical bonding, namely: metallic Na, molecular Cl 2 , and ionic rocksalt NaCl, whose electronic charge densities are reported in Figure 1. As discussed in the Introduction, the significantly different features in the electronic structure expectedly require a different support and hence a specific basis set. The geometries adopted are fully reported in the Supporting Information and have been obtained from experimental references in the literature. 47−49 In Table 3 we see that for the Cl 2 molecular crystal, not unexpectedly, the original def2 basis set undergoes very little modifications when optimized in the solid. Actually, it performs much better than the pob-TZVP basis (see Table  5) with the total energy being 0.1 Ha lower. The removal of the outermost p-function in the pob basis sets leads to an overall decrease of the exponents of the remaining functions that partly compensates the contribution to the total energy of the missing function. If one includes the outermost p-function from the dcm basis set, a further energy lowering of 13 and 17 mE h is observed for the pob and pob-rev2 basis set, respectively. However, this is not enough to reach the final energy of solid Cl 2 as obtained with the optimized dcm basis set thus showing the crucial role of the outermost p-function.
The dcm[NaCl] basis for Cl, optimized in the rocksalt structure, features significantly more contracted exponents, as far as s-and p-functions are concerned, while the d exponent becomes more diffuse. As reported in Table 4, a stronger contraction is observed in exponents of the s-type orbitals in going from the molecular def2 to the bulk metal and then the ionic NaCl. In this case we had to remove the most diffuse pfunction (0.03 au) in order to ensure convergence, but at difference with the pob-TZVP case, we were able to keep all the d-functions in.
As shown in Table 5 it can be seen that in all cases dcmenergies are significantly lower than pob-ones, and quite surprisingly the dcm[Cl 2 ]-TZVP and dcm[Na]-TZVP basis sets seem to perform well also in the ionic case.
Such basis sets effects are also reflected in geometry optimizations. In Table 6 we report the optimized lattice parameters obtained with the different basis sets. It is seen that the dcm-basis sets lead in all cases to an expanded volume with respect to the pob-ones and in the case of Na and NaCl also to a better agreement with experiment at the PBE0 level. In the molecular crystal Cl 2 , dispersion effects act as a key role, hence the plain PBE0 leads to an excessively expanded volume when the dcm[Cl 2 ] basis is used, while the introduction of -D3 dispersion correction restores a more correct description. It is reasonable to assume that the volume expansion associated with the dcm[Cl 2 ] basis is related to a mitigation of BSSE effects−which usually act as spurious dispersion.
3.3. Use of Large, Extended Basis Sets. Solid LiH is a rather standard benchmark for methods assessment in the solid state. Recently 24,26 lithium hydride has been used as a benchmark for estimating the Hartree−Fock basis set limit compared with results from different approaches. 23,50,51 The case of LiH, similarly as NaCl, poses certain difficulties since standard molecular basis sets are designed for neutral atoms, not ions, hence inapplicable to bulk ionic crystals without modification. We optimized the basis set series def2-SVP/def2-TZVP/def2-QZVP with our BDIIS algorithm, obtaining the    Figure 4 we compare such energies with previous data from the literature also obtained with the CRYSTAL code. It is seen that with the quadruple basis we reach a value that is very close to that of ref 26 (i.e., −8.06475 E h ), where a much larger basis set was used. This last result was already close to the CBS limit compared to methods employing different basis set types. 26 If the dcm-TZVP and dcm-QZVP total energies are used to estimate the HF complete basis set (CBS) limit by using a twopoint extrapolation scheme based on an exponential formula, a value of −8.065089 is attained. Notably, this energy limit is even lower than the one reached by Usvyat and co-workers 26 by 0.3 mE h . When using the CBS energy for the atoms, 52 the cohesive energy is then −3.60 eV in very good agreement with results from different theoretical approaches. 23,50,51 Similarly, we have optimized a quadruple-ζ basis for diamond and graphene. The original def2-QZVP basis does not allow convergence in either case, while with the basis sets as reported in Table 7 the energies of −76.165178 and −76.174386 E h are obtained for the two systems at the PBE level. The latter value we believe to be close to basis set completeness. Extrapolation to the CBS limit leads to a value of −76.167396 and −76.177298 E h , respectively.
In Table 7 we report the reoptimized exponents with respect to def2-QZVP basis sets−all other functions are the same as in the molecular basis set. It is worth noting that g-type functions were also included in the basis set as they were recently made available in the development version of the CRYSTAL code. 53 BSSE effects are reduced much more considerably by the increasing of the basis set quality, rather than by the optimization of the exponents, so that BSSE is quite similar for pob-or dcm-basis sets.
For diamond, we have also calculated the Hartree−Fock CBS limit by using the dcm-TZVP and dcm-QZVP basis sets.     54 Interestingly, results show that the basis sets optimized with the PBE functional can also be used for HF even if a tighter setting of the computational parameters is required (see the Supporting Information).
In Figure 5 we compare, for graphene, the electronic band structure computed with Gaussian basis sets with the bands for the same systems as obtained from a plane wave code 55 using a considerably high cutoff. It is evident that the bands at the triple-ζ level are different from the reference ones, specially in the Γ and M points of the Brillouin zone. Nevertheless, the dcm[C graph ]-TZVP performs better than the pob-TZVP. A considerably better agreement is attained by using the dcm[C graph ]-QZVP basis (right panel of Figure 5). We believe this is strong evidence of the possibility of reaching converged results with Gaussian basis sets and the effectiveness of a system-specific optimization scheme.

CONCLUSIONS
In the present work, we have developed a basis set optimizer based on the DIIS algorithm that minimizes the total energy of the system constrained to keep the condition number of the overlap matrix as small as possible in a similar approach as proposed by VandeVondele et al. 40 The latter constraint acts as a pivot in the optimization of the basis set and prevents the lowest exponents of the basis set to decrease too much thus reducing the risk of linear dependency and numerical instability. This is particularly important in solid-state calculations where the use of atom-centered diffuse functions is more delicate and sometime useless.
We have then shown that the proposed method is quite effective for solid-state calculations and allows for an easy optimization of basis sets not only of triple-ζ quality but even of quadruple-ζ size. Furthermore, we have demonstrated that the BDIIS method can be used to obtain basis sets for solids of consistent quality as molecules without pruning the original basis sets. Results for simple solids as diamond and graphene for which the definition of an appropriate and systemconsistent basis set is uglily difficult are very promising. Also, the possibility of employing basis sets specifically calibrated on a given system allowed us to easily reach the HF complete basis set limit for LiH which has been a long debated issue and for diamond.
While reasonable questions can be raised about the transferability of such optimized basis sets from one method to another, we have seen that a basis set optimized, say, with PBE is very close to convergence when inserted in HF or PBE0. For our diamond test case the energy with such basis was only a few μE h away from the minimum when transferred from one method to another.
The evidence of the excellent performance of the BDIIS method paves the way for a careful definition of system-specific basis sets, as a viable alternative to all-purpose basis sets. Nevetheless, it could be employed for a more extensive work that would permit the creation of all-purpose basis set families for a larger set of atomic species. Furthermore, the algorithm here described could be very useful to optimize basis sets for post-HF correlation methods 20,56 as well as for response properties. 57 The authors declare no competing financial interest. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article ■ ACKNOWLEDGMENTS