Hybrid grid/basis set discretizations of the Schrödinger equation

Most numerical simulations of the Schrödinger cannot work directly in the continuum; they require a discretization of some sort that maps the continuum to a finite number of degrees of freedom. For example, in the many thousands of electronic structure calculations performed annually for molecules and solids, basis sets are usually used. For molecules, these are usually atom centered functions composed of linear combinations of Gaussian type functions, for which analytic integrals are used to quickly generate the discrete Hamiltonian terms. For solids, plane wave basis sets are often used, which also have analytic integrals. These basis set methods have been developed and refined over most of a century, and the software to use them has been improved to a remarkable degree. Nevertheless, a variety of alternative discretization methods such as wavelet bases,^1–3 adapted grids,^4,5 finite elements,⁶ and sinc-function bases⁷ have continued to be developed. These alternatives have been pursued to address real weaknesses in existing discretizations, but Gaussian and plane wave bases are still used for most calculations.

What is the motivation to try to make an improved discretization approach again? Our primary reason is that some of the most promising new strong-correlation methods for electronic structures—the density matrix renormalization group (DMRG)^8–10 and related tensor network methods¹¹—depend strongly on the discretization. These methods have strong requirements for locality in the basis, due to the low-entanglement nature of their approximations. One can use localized Gaussian basis sets for DMRG, but this combination does not perform nearly as well as DMRG for simple lattice models, like the Hubbard model. The ideal basis for a DMRG or tensor network electronic structure calculation would be some sort of hypothetical real space grid with a very modest number of grid points, comparable to the number of functions in an atom-centered Gaussian basis, and giving comparable accuracy. Besides the DMRG, some other approaches in the electronic structure, such as linear scaling methods, also benefit or require high locality and grid or grid-like representations, and this has stimulated substantial work in this general area. However, the specific requirements of DMRG and related methods, which differ from other correlation methods, leave room for new approaches distinct from previous work.

Our approach has two key elements. The first is a novel type of basis function, called a gausslet, related to the scaling functions of wavelet approaches, designed to live on a uniform grid. Several features of these functions are particularly nice for electronic structures, particularly that they are orthogonal, localized, smooth, symmetric, and composed of a sum of Gaussians. The second element is a set of diagonal approximations for the Hamiltonian elements. These sort of diagonal approximations have previously been developed for sinc-function bases,⁷ but we show that these approximations also arise in another context, in terms of wavelet-type bases with certain moment properties, such as our gausslets. The gausslets are much more localized than sinc functions, making them better suited for DMRG. In making a diagonal approximation for the two-electron interaction, the number of terms (integrals) needed drops from N⁴ to N², where N is the number of basis functions. This reduction is not a large-N asymptotic feature—it appears for any N with no large coefficient in front of N². There are other recent approaches, for example, utilizing density fitting and tensor factorization,^12,13 that also produce O(N²) representations of the two-electron interaction. The DMRG might be adapted to use these representations; our approach would not require any changes to the standard DMRG. The gausslets and diagonal approximations combine to take us a substantial way towards the performance of the ideal very coarse but accurate real space grid.

Section II gives a more detailed overview of the problem and the approach. Section III discusses grids of Gaussians as bases, which have nearly ideal completeness properties but suffer from near linear dependence. Section IV discusses ternary wavelet transformations (WTs) that are used in Sec. V to derive gausslets. Section VI presents diagonal approximations, and Sec. VII has conclusions and a discussion of further directions. Throughout we give numerical examples in 1D.

Consider an arbitrary basis set approach. The Hamiltonian is an operator with the coefficients of the operator terms defined by integrals over the basis. By far, the most complex part of the Hamiltonian is the two-electron Coulomb interaction operator, parameterized by a four-dimensional array or tensor

where ϕ_i(r) make up the set of N basis functions, i = 1, …, N. In some simpler approximations (e.g., the local density approximation of density functional theory), one can bypass the use of V_ijkl in favor of, say, solving the Poisson equation, but usually for more accurate methods treating electron-electron correlation explicitly, one cannot avoid dealing with V_ijkl, which has N⁴ elements.

In contrast, a simple cubic grid discretization gives a representation of the electron-electron interaction which scales only as N², where now N is the number of grid points. In such an approach, one could use finite-difference approximations to represent kinetic energy derivatives. The nucleus-electron and electron-electron interactions would be evaluated point-wise, with the two-electron Hamiltonian terms taking the form $Vijn^in^j$⁠, where $n^i$ is the density operator on site i, and $Vij=|r→i−r→j|−1$ (with a suitable alteration at i = j).

Simple grids like this are rarely used for electronic structures because N for an accurate 3D grid would be much bigger than for a typical basis. Consider a single atom: the grid spacing needs to be set to resolve the behavior near the nucleus, resulting in many grid points to describe the tails of the wavefunction far from the nucleus. In contrast, in a Gaussian basis, the tails of the wavefunction can be described using just a few basis functions. For example, in an early effort by the author using a uniform 3D grid, utilizing finite elements rather than finite differences, as many as 10⁵ grid points were needed for accurate results for the molecule.⁶ Much more efficient approaches use adapted grids, putting more points near the nuclei.⁴ However, it is difficult to make the number of grid points nearly as small as the number of basis functions in a Gaussian-type basis. Wavelet bases are another approach with much of the locality of a grid and with adaptable increased resolution near nuclei, but again the number of functions tends to be substantially larger than that in Gaussian bases. The larger number of functions is a particular problem for methods which treat correlation especially accurately.

Grid-like methods have a fundamental advantage for use with DMRG and other recently developed tensor network methods: these methods are based on the low entanglement of ground states when expressed in a local real-space basis. The entanglement of ground states is governed by the area law,^14,15 which is specific to localized real space bases. In a delocalized basis, a volume law of entanglement holds instead, requiring the number of states m kept to grow exponentially for fixed accuracy. Thus the DMRG greatly prefers a real-space local representation. In addition, the calculation time for the DMRG depends strongly on the number of two-electron interaction terms, so the N⁴ scaling of V_ijkl is a major drawback. The standard DMRG approach for molecules in a Gaussian basis utilizes a standard basis localization method before the DMRG is used. This localization is less than ideal, and except in certain treatments of long chains, all N⁴ two-electron terms are kept.

The very recently developed sliced basis DMRG (SBDMRG)¹⁶ approach uses a finite-difference grid in one direction (the long direction in a chain) and 2D Gaussians for the transverse directions. For chain systems, the SBDMRG grid gives a reduction to O(N²) two-electron terms, where N scales linear with the length. SBDMRG also utilizes a compression algorithm for the long range terms, yielding an O(N) calculation time, and chains of up to 1000 hydrogen atoms have been studied in the strongly correlated stretched-bond regime. While arbitrary molecules can, in principle, be treated with SBDMRG, it is particularly suited for long chains, even more than the standard quantum chemistry DMRG. One motivation for this work is to find a good way to increase the grid spacing in the SBDMRG without loss of accuracy, for faster and more efficient SBDMRG calculations. However, a more important motivation is to take advantages of locality in all three dimensions.

Here we present techniques that do take us closer to an ideal combination of grid and basis discretizations. Our approach is related to orthogonal wavelet bases. Wavelets have and continue to be used successfully in electronic structure calculations, but our goals are different from many existing uses. A conventional wavelet basis can be used in a way that has precise control over the accuracy, with systematically improvable accuracy, in, say, a density functional theory calculation.¹ These approaches make no use of standard Gaussian bases which give an excellent description of core electrons with a very small number of functions. Thus the price of the precise control of accuracy is a basis that may have more than 1000 functions for a small molecule, which is not practical for most beyond-HF correlation calculations.

Our interest is in these correlation calculations, where one is interested in accuracies only up to about 1 mH and for which Gaussian basis set descriptions of cores are fine, at least for energy differences. Thus we consider approaches where we use wavelet techniques in a very restricted way, combining them with standard Gaussian basis sets. A key consideration in developing such an approach is the calculation of integrals. To make this easy, we construct functions similar to the scaling functions of wavelet transforms, but out of an equal-spaced array of identical-width Gaussians. The resulting functions called gausslets have a number of significant advantages for electronic structure calculations.

Gausslets live on a uniform grid, with one function per grid point. They are defined in 1D, but in 2 or 3D, one simply takes products of the 1D functions. Gausslets are orthogonal and symmetric and almost compact in both real and momentum space, in the same sense that Gaussians are. They have excellent completeness properties for representing smooth functions. All standard integrals have simple analytic forms although one has to contract over the underlying Gaussians. A gausslet is shown in Fig. 1, together with seven of the Gaussians of which it is composed. A 1D gausslet basis would be composed of this function plus integer translations, and the whole basis can be scaled to a grid spacing a.

Gausslets are constructed using a recently developed class of wavelet transforms which use a factor of three-scaling transformation rather than the usual factor of two, in order to make them exactly symmetric. They have an additional very important property: they integrate like a δ-function, for any low-order polynomial p(x),

This, in turn translates to reducing the two-electron integrals from N⁴ to N², V_ijkl → $Ṽij$⁠, provided the interaction is smooth. One way to make the interaction smooth would be to replace the Coulomb electron–electron interaction with a two-electron pseudopotential.¹⁷ 

A motivation for the gausslet development was to make a basis that was easy to use in combination with standard Gaussian bases. In particular, if one thinks of a standard grid as being defined by finite differences, it is not at all clear how one can add, say, an extra atom-centered narrow-width Gaussian to the representation to describe a nuclear cusp. But as soon as one chooses to attach basis functions to each grid point, then it becomes clear how one can add additional non-grid functions. This would be true also for a grid of sinc functions, and with FFT methods one can do integrals with sinc functions quickly. However, we also wish to orthogonalize the combined basis (where either gausslets or sinc functions are orthogonal among themselves). Ideally, we would also want to be able to orthogonalize separately around each atom. This gives very localized bases, such as gausslets, an advantage. The number of gausslet-Gaussian integrals is reduced by the locality, and their evaluation is made easier by the fact that the gausslets are composed of Gaussians. These practical issues involving combining Gaussians and gausslets, and doing the needed integrals, are important but are beyond the scope of the current paper.

The first step in constructing gausslets is understanding the properties of arrays of equally spaced, equal-width Gaussians.

Arrays of Gaussians with identical widths are particularly convenient to use in constructing basis functions because of their analytic integrals, their smoothness, their completeness, and also because the product space is greatly reduced and convenient. The product of a Gaussian of width $w$ located at i, and another of width $w$ at j, is a Gaussian centered at (i + j)/2 of width $w/2$⁠; the set of all products of a grid of N such Gaussians is a half-spaced grid with roughly 2N functions, rather than N² functions. The product space plays a central role in defining Hamiltonian matrix elements, so this simplification can significantly improve computational efficiency. (Polynomials also have similar features of their product space, namely, that a product of two-order N polynomials is an order 2N polynomial. This representation is used in the exact tensor hypercontraction approach.¹³)

The completeness properties of arrays of Gaussians are interesting. As is common in numerical analysis, we define completeness in terms of representing low order polynomials. A grid of basis functions with good polynomial completeness is excellent at representing arbitrary smooth functions. Consider a unit-spaced 1D grid, and on each integer point j, put a Gaussian with width $w$⁠,

To see the completeness of the set {g_j}, consider the sum, shown in Fig. 2(a),

Since g_j are identical except for translation, the representation within this basis of a constant function must be proportional to C(x). If C(x) is not nearly constant, the basis is very poor—it cannot even represent a constant.

We will show that C(x) is very nearly constant for sufficiently large $w$⁠. Note that C is periodic, C(x − 1) = C(x), with period 1. Expand it in a Fourier series, with coefficients c(k = 2πn), where n is an integer

The largest deviation from constancy comes from the smallest nonzero k in the series, namely, k = ±2π, for which the factor $exp(−12k2w2)$ is exp(−2π²w²). For $w=2$⁠, this factor is less than 10⁻¹⁷. For $w$ = 1, it is less than 10⁻⁸. Thus for modest width $w$⁠, C(x) is exactly constant to double precision accuracy (15 or 16 digits). In what follows, we choose $w$ = 1 since this lack of completeness at the 10⁻⁸ level appears in the wavefunction. According to standard variational arguments, it translates to errors in the energy of order 10⁻¹⁶. Note that higher terms in the series have much smaller coefficients so that not only C(x) is nearly constant but also a number of higher derivatives are nearly zero. We will not try to make this particularly quantitative since we can test any such property with a trivial numerical evaluation.

The flatness of C(x)—for the special case of Gaussians—also implies higher order polynomial completeness. In particular, we will show that the near constancy of C(x) implies that, for n ≥ 1, ∑_j jⁿg_j(x) is very nearly a polynomial in x with order n. (We have already shown it is true for n = 0.) If some linear combination of g_j can represent an order-m polynomial for 0 ≤ m ≤ n, then one can find a suitable linear combination that can represent any polynomial up to order n. Gaussians have the property that their derivatives, to all orders, can be written as the same Gaussian times a polynomial. In particular, let P(x, n) denote a polynomial in x of degree n; the nth derivative of g_j(x) can be written as

where C_l,m = 0 if l + m > n. The nth derivative of C(x) is

In order to prove our result by induction, we pick out the l = n term, for which m = 0, and assume that for l < n, ∑_jj^lg_j(x) is a polynomial P_l(x). Then

Since C_l,m = 0 if l + m > n, this means that if ∑_jj^lg_j(x) is very nearly an order-l polynomial for all l < n, then it is an order n polynomial for l = n. By induction, this is true for any n. This result breaks down when C⁽ⁿ⁾(x) stops being zero for large n. For $w$ ≳ 1, the breakdown does not occur very quickly, and an array of Gaussians has approximate polynomial completeness to high order.

The weakness of a grid of Gaussian as a basis is its lack of orthogonality, and that orthogonalizing the functions involves a fairly singular matrix, associated with a near lack of linear independence. Let S( j, k) = S( j − k) = 〈g_j|g_k〉 be the overlap matrix of {g_j}. We can form an orthonormal set of functions (all related by integer translations) using q(j − k) ≡ S^−1/2(j, k) as

However, S becomes increasingly singular as $w$ increases. The eigenvectors of S are plane waves, and the most singular point is at momentum π; this corresponds to near cancellation of g_j. For example, for $w=2$⁠, the unnormalized fit to a momentum π plane wave, ∑_j(−1)^jg_j(x), is only of order 10⁻⁴ in magnitude. The near singularity of S means that G_j have long tails and widely varying coefficients for the component Gaussians, making them computationally poorly behaved.

Fortunately the orthogonalization does not spoil the completeness nor would a generic convolution of the form Eq. (8), transforming g → G using an arbitrary vector q. Consider

The term in square brackets is a polynomial in k of degree n, so the entire right side is a polynomial in x of degree n. By the same reasoning as above, this implies that G_j have the same excellent completeness properties as g_j.

In summary, a grid of Gaussians has excellent completeness, as well as several other very desirable features, but poor orthogonality and linear independence properties. Fortunately, as we show in Secs. IV and V, these faults can be fixed using wavelet technology.

In conventional compact orthogonal wavelet theory, as pioneered largely by Daubechies,¹⁸ it is not possible to have a symmetric, compact, orthogonal wavelet transformation (WT). However, with ternary WTs, which change scales by a factor of 3 instead of 2, symmetry, orthogonality, and compactness are compatible. Recently Evenbly and White (E and W)¹⁹ introduced new families of ternary WTs with excellent completeness, compactness, and smoothness properties, which we will make use of and extend. A very carefully chosen set of coefficients c_k, − m ≤ k ≤ m, with c_−k = c_k, defines a symmetric wavelet transform. Given a single function f(x), we define all integer translations f_j(x) = f(x − j) to form a basis set that lives on an integer grid. Then the wavelet transform produces a new basis set, also living on a grid with unit spacing, defined by a function f′, defined by

The scaling function of the WT is the fixed point of Eq. (10). The fixed point exists provided the initial f satisfies some simple conditions, most importantly that they sum to a constant (exactly), ∑_jf_j(x) = const. Thus any functions defined by Gaussians with finite $w$ are technically not suitable for generating the fixed point. However, if one only wants to apply the WT a modest number of times, Gaussian derived functions can be an excellent starting point. Since we do not plan to use the wavelets to represent sharp core functions (instead using standard Gaussians), a full multiscale resolution analysis is not needed, and for our uses, the fixed point is entirely unnecessary.

E and W showed that WTs have a direct correspondence with quantum circuit theory. The circuit corresponding to a WT is defined by a small number of angles θ_k; each angle defines a unitary “gate,” which is a 2 × 2 or 3 × 3 unitary matrix. The unitarity of the circuit is independent of the values of θ_k, removing an annoying set of constraints when optimizing a WT. The circuits express symmetry much more naturally than in conventional approaches. Some of the symmetric ternary WTs constructed by E and W appear to have better properties than any previous such WTs.

Besides the transformation of the scaling function in Eq. (10), the WT produces two wavelet functions per scaling function. The scaling function captures the lowest momentum behaviour, while the wavelets describe high momentum. These wavelets are less central to our discussion here.

In the terminology of E and W, we consider here only type I site-centered symmetric ternary wavelets. These are characterized by the number of low and high frequency moments associated with the scaling functions. We make this more specific here with the notation W_nlh to describe the WT with n angles, l low moments, and h high moments. The number n also gives the number of layers of the circuit, fixing the range of c_k, specifically from −(3n − 2) to 3n − 2, with c_−k = c_k. The number l gives the completeness of the scaling functions; more specifically, l − 1 is the maximum polynomial order the scaling functions fit exactly, which is imposed by making the wavelets orthogonal to any lth degree polynomial. Similarly, h controls the smoothness; specifically, h − 1 gives the number of sign-flipped polynomials (−1)^kk^h that the scaling functions are orthogonal to; it also means that, the Fourier transform of c_k has vanishing value and h − 1 vanishing derivatives at maximum momentum π.

A related and important property of the corresponding scaling functions is that they integrate polynomials like a δ-function. This is associated with a property of c_k,

for m = 0, …, p − 1 for some p. (The evenness of c_k makes this statement trivial for odd m.) The δ-function order p is related to the other moments but not in a simple way: we have found WTs with the same n, l, and h, but with different p. Usually, however, a large p appears “for free” from optimizing l and h.

E and W give just a few examples of these types of WTs and not to full double precision. Determination of W_nlh involves a nontrivial nonlinear optimization; we have found W_nlh to high precision both for the examples of E and W and a number of additional WTs with higher order. All useful WTs had even n, and it is difficult to converge useful WT for higher n. Angles for the most useful WTs are listed in Table I. These emphasize completeness over smoothness, but not completely, mostly setting h = 2 (which is only one constraint since the order-1 high moment is automatically satisfied because of symmetry). The scaling functions with h = 2 are nicely smooth; for h = 0, they are much more irregular. This leaves n − 1 degrees of freedom for completeness, i.e., l = n − 1. Below we construct gausslets out of W₄₃₂, W₆₅₂, W₈₇₂, and W₁₀₉₂. These WTs all have impressively high p’s given by a simple formula: p = 2l.

A WT can be applied to an array of Gaussians, producing a new basis where the functions are formed from sums of Gaussians. The completeness properties of the Gaussians are transferred by the WT to the new basis, up to the completeness order l. To see this, let i index a grid with spacing 1, and j (an integer) run over a grid with spacing 1/3. Let f_i(j) = c_j−3i; we can think of f_i as being a discrete basis function representing functions living on the 1/3 grid. (Usually in wavelet theory, one thinks of an infinite sequence of ever smaller grids, which represents the continuum in the limit of infinitesimal grid spacing. In contrast, we apply the WT to a grid of functions which are already continuous functions of x.) The WTs have a discrete completeness property that means we can write

for m < l. Multiplying by g_j(x) on both sides and summing over j, we find that linear combinations

can represent arbitrary low order polynomials. Again, this implies that suitable linear combinations of s_i(x) can represent any low order polynomial. However, they are not orthogonal, which we address in Sec. V.

We wish to use these properties of WTs and Gaussian grids to define convenient, orthogonal, complete functions defined as sums of Gaussians. There are several possible ways to proceed. One way starts with s_i(x). s_i(x) are not orthogonal, but their overlap matrix S_s is much less singular than that of g_j(x). s_i(x) look like somewhat broadened wavelet scaling functions, with oscillations that partially cancel off-diagonal elements of the overlap matrix. (Further applications of a WT would take them closer to orthogonality.) One can use $Ss−1/2$ to orthogonalize s_i. The resulting functions are close to what we want, but $Ss−1/2$ introduces moderately long tails in the final functions. One can then apply a WT to these orthogonal functions, which would substantially decrease the tails, and use those functions. The main drawback is that the underlying Gaussian grid would have a spacing of 1/9 that of the final basis functions. The more the Gaussians making up a basis function, the more is the computational work to use them, so we have developed a way to maintain the minimal 1/3 spacing of the underlying Gaussians without the long tails.

It is possible to partially orthogonalize g_j(x) so that the application of a particular WT to the partially orthogonalized functions produces fully orthonormal functions. Let function G_j(x) be defined by Eq. (8), for some vector q(k), and then apply a WT [defined by f_i(j)],

We know that q(k) exists that makes the $Gi(x)$ orthonormal—we can use the overlap matrix S_s to make q(k), which completely orthonormalizes G_j(x), and thus $Gi(x)$⁠. In that case, q(k) falls to about 10⁻¹² near k = 110, so it is rather nonlocal, making $Gi$ have long tails. It is possible to find a much more compact q(k) if it is optimized for a particular WT, requiring orthogonality only of $Gi$ and not both G_j and $Gi$⁠. The q convolution and the WT, each take care of parts of the orthogonalization. A typical q is shown in Fig. 3. Finding such q is a rather tedious nonlinear optimization problem, with many local minima. We implemented this optimization using a Nelder Mead algorithm, using restarts and high precision to deal with local minima, with the minimizations sometimes taking several days, and only achieving approximate double-precision orthogonality. A final orthogonalization step was performed to make the functions orthogonal to very high precision, with minimal impact on locality, using the inverse square root of the overlap matrix. Fortunately, using these functions does not require that these optimizations be repeated; all one needs are the final coefficients of the Gaussians, and the properties of the functions are easily verified.

This approach produces functions—gausslets—with a ratio of 3 between the spacing of the functions and the underlying Gaussians. Normally with WTs, there is a tradeoff between order and compactness. The presence of the underlying Gaussians means that there is an inherent limit to the size of a “1/3-gausslet” which is greater than the width of the fixed point scaling function of the WT used to make the gausslet. However, the gausslets are substantially smoother. We present coefficients for gausslets $G4$⁠, $G6$⁠, $G8$⁠, and $G10$⁠, where the index signifies n of the underlying WT, W_n, which has l = n − 1 and h = 2. The gausslet $G6$ is shown in Fig. 1; the others look similar, with slightly increasing width with the order. The coefficients b_j of these gausslets in terms of the underlying Gaussians are given in the  Appendix in Tables II–V.

A simple test of the completeness of the gausslets is shown in Fig. 4. The gausslets are fitted to polynomials of various orders and the root mean square errors over the unit interval are shown. (The gausslets involved in the fitting extend well outside this interval because of the small but finite tails of the gausslets.) The fitting coefficients were approximated using the δ-function property, as T_n(i/10) for the gausslet centered at i, so this test of the fitting also tests the δ-function property. Since the order of the δ-function property is much higher than for fitting property, one sees only the errors associated with the fitting.

One may wish to have more compact gausslets, paying the price of an underlying finer array of Gaussians. Suitable 1/9, 1/27, etc., gausslets are easily derived by applying WTs to one of the $Gn$⁠. We label these gausslets as in the following example: $G864$ is formed from $G8$ and applying W₆ and then W₄. These constructions also generate useful wavelet-like functions, and a multiresolution basis can be formed for a modest range of scales, with each length scale having different functions. If one successively applies W_n, the range of the functions decreases, and the basis functions become close to the scaling function of W_n. Alternatively, one can successively decrease the order over several steps, which produces very localized functions which are still smooth with a finite size underlying Gaussian array. An example of this is shown in Fig. 5. A multiscale basis can be useful, but in many cases, a better alternative that produces smaller bases is to add atom-centered Gaussian functions to the gausslet basis. The extra functions can be orthogonalized to the gausslets and themselves. A simple 1D example of this approach is given below.

As a first test of the use of gausslets to solve the Schrödinger equation, we apply them to two different simple 1D smooth potentials; see Fig. 6. The first, the Pöschl-Teller potential, is exactly solvable, with energy E = −1/2. The second, the soft Coulomb (SC) potential, has been used as a 1D model with properties such as correlation energies roughly similar to that of real 3D molecules.^20,21 Here the SC potential would be that of pseudo-hydrogen. Each potential is centered at two positions to demonstrate that the grid points need not be aligned with the potential centers. To define the Hamiltonian matrices, we first evaluate the Hamiltonian matrix elements for the underlying Gaussian grid with spacing a/3, analytically for the kinetic energy and numerically for the potentials. Then we transform to the gausslet basis using b_j. Both of these potentials have widths of order 1, and we see that very high accuracy is easily achievable for a ∼ 0.2. More important for our goals is that one achieves accuracies of O(10⁻³) for a ∼ 1. Large a values are essential for applications to 3D.

A key advantage of a basis set approach is that one can add functions to it to represent singularities or other sharp features. As a test of this, shown in Fig. 7, we consider a potential with a smooth part plus an off-center delta function, $V(x)=−exp(−12x2)−δ(x−2.7)$⁠. Particularly in three dimensions, this sort of potential is difficult with a grid. Our basis consists of an array of gausslets plus, to represent the singularity, the solution to the delta function alone, exp(−|x − 2.7|). The exponential is orthogonalized to the gausslets to make a fully orthogonal basis. To carry out the calculations, the exponential is first represented very accurately as a sum of a few hundred Guassians, the coefficients coming from a discretization of an integral. (In a 3D calculation, to describe, say, a 1S function, one would use a sum of a few Gaussians from a standard basis set.) Then all the potential terms are sums of analytic Gaussian integrals. The hybrid Gaussian/exponential approach converges very rapidly, with errors below about 10⁻³ for a < 1.5 with $G8$ and $G10$⁠. The energies for the last two points, a = 0.1 and a = 0.2, agree to eight digits, E = −0.66 144 716. For comparison, we implemented a naive grid, without trying to put a grid point exactly on the singularity. The grid results show substantial oscillations as a function of a. The finest grid result, a = 10⁻⁴, gives E ≈ −0.6614, in agreement with the gausslet result but much less accurate.

A crucial feature of finite difference grid approximations is the simple form of the two-electron interaction, $Vijn^in^j$⁠. A similar property holds for the one electron nucleus-electron potential terms, which have the simple form $Uin^i$⁠. This latter property is not so important in itself since the extra computation time to deal with a U_ij form of the potential is minor, but it does serve another useful purpose: approximations which make U diagonal should generally do the same for V since the two-electron interaction acts on one electron with a single-particle potential determined by the locations of all the other electrons. We will consider the simpler single particle potential first, in one dimension. Diagonal representations such as these are a key property of bases made from the sinc function.⁷ Here we find similar approximations for gausslets, implying that the high nonlocality of the sinc is not needed for this sort of approximation. Although we focus on the gausslets, similar approximations would work for the ternary wavelet scaling functions and for certain types of traditional wavelet scaling functions with the δ-function property, particularly coiflets.¹⁸ Devising methods to reduce V_ijkl to matrices or tensors of size O(N²), using tools such as density fitting and tensor factorization, is also a central focus of other modern approaches.^12,13

Let U(x) act on a continuum wavefunction ψ(x) to give another wavefunction ϕ(x): ϕ(x) = U(x)ψ(x). We are interested in the coefficient of ϕ for basis function $Gi$⁠, written using the expansion of ψ in terms of $Gj$⁠,

This is the conventional non-diagonal form, with the last integral defining U_ij. Now assume that each $Gi$ integrates like a weighted delta function, with location x_i,

for any smooth function f(x), with

Then

Thus we have a diagonal “point” approximation for U,

Note that the locality of the gausslets means that U_ij is small for i far from j, but the diagonal approximation does not rely on this. In particular, we disregard the near-neighbor j = i ± 1 terms even though they are not small.

There are two other closely related diagonal approximations. First, if the δ-function property is only approximate, we may prefer to replace U(x_i) with its overlap with $Gi$⁠,

One might hope that this “integral” approximation could be more accurate since it averages the potential over a finite range, but one needs tests to see if it actually is an improvement. Second, let us assume that ${Gi}$ can exactly represent a constant function over the range of interest. Then $wk$ is then the expansion coefficient for $Gk$ to represent the identity function, and

Inserting this into Eq. (20) gives the “summed” approximation

For the special case where $Gi$ are a uniform grid of gausslets, $wk$ = $wi$⁠, and as the grid spacing goes to zero, smooth functions have ψ_k ≈ ψ_i in the close vicinity of a point i. In this case, ∑_kU_ikψ_k ≈ ψ_i∑_kU_ik, providing another rough justification for Eq. (22). (This idea was the original inspiration for these diagonal approximations.)

Tests of the diagonal approximations are shown for the soft Coulomb potential in Fig. 8. One sees that all of the diagonal approximations become very precise for small a. In fact, although one cannot see it in this figure, the diagonal approximation curves and the non-diagonal full matrix curve come together faster as a function of a than any of the curves converge to the exact result. This is because the δ-function property parameter p is twice that of the completeness parameter l. This means that at small a, there is no point in using the full matrix—the result might not be exact, but the diagonal approximation does not contribute significantly to any error. At larger a, the diagonal approximations introduce significant errors. The integral approximation (or equivalently for this case, the summed approximation) is significantly better than the point approximation.

These three approximations translate immediately into two-electron approximations for V. The pointwise evaluation is

The integral approximation replaces V (x_i, x_j) in Eq. (23) with

The summed approximation replaces V(x_i, x_j) with

Because of the big difference in computational scaling associated with two-electron versus one-electron terms, it is sensible to only make the diagonal approximation for the two-electron terms, using the full matrices for the one electron terms.

As a simple interacting example calculation, we perform an exact diagonalization for a 1D helium atom in a soft Coulomb potential, with a basis of gausslets ranging from −L to L. The results are shown in Fig. 9. According to Ref. 21, the exact energy from a grid DMRG calculation is −2.238. The two full/integral curves, which completely overlap, demonstrate that L = 7 gives excellent convergence in L. All the diagonal approximations perform excellently at smaller a. For larger a, they still behave well; the full/integral approximation stays within 1 mH almost up to a = 1. The point diagonal approximations are somewhat less accurate at larger a, but the point/point approximation is especially convenient for small a since it is as accurate as the other approximations and requires no integral evaluations at all. For this case, we consider a = 0.1 and L = 15, which gave 601 basis functions, to push for very high accuracy; we obtain E = −2.238 257 824, which we believe is correct for all the digits shown. The setup and diagonalization took less than 5 min on a desktop. With DMRG, it would be straightforward to extend these calculations to long chains.

One can consider making the diagonal approximations after first making a mean-field-like reduction of the operator. For the one electron case derived below, this gives an energy correction but does not change the Hamiltonian. For the two-electron terms, this produces altered one electron terms which change the diagonalization and might improve the results. For the one electron case, assume we have some estimate for $⟨ĉi†ĉj⟩$⁠, and let

If we now apply a diagonal approximation U_ij → $Ũii$ on the first term, we have

Suppose we find the ground state using $Ũ$ and from this compute $⟨ĉi†ĉj⟩$ and $⟨n^i⟩$⁠, and then plug this into Eq. (28). Then adding in the correction term in the parentheses in Eq. (28) is the same as evaluating the energy in the full Hamiltonian of the eigenstate from the approximate Hamiltonian. This corrected result is shown in Fig. 8. We see that the corrected result is significantly better than the uncorrected result at large a. We leave exploring these mean-field approximations for the two-electron term for future work.

A basis permitting a good diagonal approximation is special, and an orthogonal transformation of the basis will, in general, spoil the approximation. A truncation of the basis may not harm the approximation; in particular, a reduction of a particular set of orbitals into one orbital which is a linear combination from the set still allows a diagonal approximation. This is because a basis rotation takes two-particle fermion number-conserving operators (i.e., $ĉi†ĉj$ or $n^i$⁠) into similar two-particle operators, and a single orbital only has one such operator, $n^i$⁠. This applies also to the two-electron terms. If one knows the single particle reduced density matrix (RDM, $⟨ĉi†ĉj⟩$⁠), one can determine if reducing a particular set of functions to one function is a good approximation: it corresponds to the diagonal block over the sites of the RDM having just one significant eigenvalue. If so, that eigenvector would be the correct linear combination. One may expect that the outer tail regions of a molecule would satisfy this criterion since far from the nuclei, the wavefunction decays in a simple way. This suggests that one should combine sites in the tail regions into tail functions, which could lead to a large reduction in the number of basis functions. If one transforms a set of basis function into a few functions, this breaks the diagonal approximation, but it only generates a limited number of off-diagonal terms, which may be acceptable. Specifically, if basis function i is in the reduced set and p, q, and r are outside the set, terms such as $ĉi†ĉp†cqcr$ with only one (or three) operators in the set are not generated.

A very important question is how to make a diagonal or near-diagonal approximation with a basis composed of an array of gausslets with some additional non-gausslet functions, particularly to represent the sharp behaviour near nuclei. We leave this question for future studies.

Application of these ideas to sliced basis DMRG would require little development beyond what is presented here. The initial applications of SBDMRG used a finite difference grid which was adapted with a low-frequency filtering procedure to increase convergence with the grid spacing. Usually a grid spacing of 0.1 bohr was used for hydrogen chains. The grid plus filtering procedure could be replaced by a gausslet basis, with a diagonal approximation introduced for the two-electron interaction. Probably this would allow substantially higher lattice spacing, and it would smooth the way for including larger Z atoms.

More generally, for realistic 3D systems, one can form a 3D gausslet basis by making products of 1D gausslets, namely,

If one is not worried about making diagonal approximations, then one can combine gausslets with 3D Gaussians from a standard Gaussian basis, orthogonalizing them to the gausslets. The fact that gausslets are made of Gaussians would facilitate this. Note that for a uniform grid of gausslets, all the two-electron integrals between gausslets need be done only once and would be applicable to all systems, with a simple rescaling for different lattice spacings.

Diagonal approximation can be extremely effective in speeding up calculations, and so one would want to try to find approaches that allow them. The most straightforward approach would involve using a pseudopotential both of the usual one-electron type and also a two-electron pseudopotential.¹⁷ This one-electron pseudopotential means that a gausslet 3D grid, without supplementary Gaussians, would be adequate, and the two-electron pseudopotential would allow any of our diagonal approximation for the two-electron interaction. The two-electron diagonal approximation would greatly speed up DMRG and also allow other tensor network methods, such as projected entangled pair states (PEPS), to be tried. The diagonal interaction could be compressed for use in DMRG, just as is done in SBDMRG. The diagonal approximation may also speed up other correlation methods. To the best of our knowledge, the combination of both types of pseudopotentials has not been used before, but it seems reasonably straightforward. The key question would be: what grid spacing is needed for reasonable accuracy?

For all electron calculations combining gausslets and 3D Gaussians, two-electron diagonal approximations would need to be tested and developed since the extra functions do not fit within the δ-function framework we used to derive diagonal approximations. Nevertheless, one might find approximations with acceptable accuracy. One could also consider whether diagonal approximations could be useful even in a standard Gaussian basis which has been localized by a standard method. Recently, Baker, Burke, and White²² have proposed “wavelet localization” (WL), where an auxiliary wavelet basis is used to localize an existing delocalized basis, at the cost of a modest increase in the number of functions in the basis. One could equally use gausslets to perform wavelet localization. One scheme would be to WL to localize standard Gaussian atom-centered basis functions on each atom. The result would be functions which look like standard Gaussians close to their nucleus but midway to the next atom they would rapidly die off with oscillations to make them orthogonal to all functions on neighboring atoms. This might give rise to a modified diagonal approximation in which V_ijkl is nonzero only if i and l are on the same atom and j and k are on the same atom.

In conclusion, we have introduced gausslets, a new type of basis function, which combine the efficiencies of working with Gaussians, but with systematic completeness and orthogonality, while maintaining locality, symmetry, and smoothness. We have introduced diagonal approximations, which are tied to gausslets, which dramatically improve the scaling of electronic structure calculations, making them act more like a grid than a basis. Although the various tests we have performed here are in 1D, these basis functions were developed with 3D applications in mind, and we anticipate rapid development of 3D uses.

We thank Miles Stoudenmire, Glen Evenbly, Kieron Burke, Takeshi Yanai, Tom Baker, and Garnet Chan for helpful conversations. We acknowledge support from the Simons Foundation through the Many-Electron Collaboration and from the U.S. Department of Energy, Office of Science, and Basic Energy Sciences under Award No. DE-SC008696.

Here we give coefficients for the wavelets and gausslets described here. First, we present the coefficients of the gausslets, which are defined by

b_j for gausslets with even order 4 through 10 are given in Tables II–V.

In Tables VI–XI, the coefficients defining ternary wavelet transforms are given. These can be used to define sequence of gausslets based on the primary gausslets given in Tables II–V. For example, applying a WT to a primary gausslet gives a gausslet with Gaussian spacing of 1/9. Combining Eqs. (A1) and (10), we have

where b_j come from the primary gausslet and c_k are the coefficients associated with the additional WT. To define wavelet-like functions coming from this transformation, we use the same formula, but with the mid or high wavelet coefficients replacing c_k. One can repeat these WTs to construct a multi-level wavelet-like basis, where the functions at each scale, in addition to different scalings, are slightly different. All the basis functions in the multi-level basis would be written in terms of a single grid of Gaussians, which would have a spacing a factor of three smaller than the finest level of wavelet-like functions.