Basis set convergence of Wilson basis functions for electronic structure

TL;DR

This paper investigates the convergence properties of Wilson basis functions, a phase-space localized basis set, for electronic structure calculations, comparing their performance with wavelet transformed Gaussians and discussing implementation strategies.

Contribution

First comprehensive study and numerical implementation of Wilson basis functions for electronic structure, including convergence analysis and comparison with gausslets.

Findings

Wilson basis functions show promising convergence properties.

Augmentation with Slater-type orbitals improves description of diffuse orbitals.

Wilson basis compares favorably with wavelet transformed Gaussians.

Abstract

There are many ways to numerically represent of chemical systems in order to compute their electronic structure. Basis functions may be localized in real-space (atomic orbitals), in momentum-space (plane waves), or in both components of phase-space. Such phase-space localized basis functions in the form of wavelets, have been used for many years in electronic structure. In this paper, we turn to a phase-space localized basis set first introduced by K. G. Wilson. We provide the first full study of this basis and its numerical implementation. To calculate electronic energies of a variety of small molecules and states, we utilize the sum-of-products form, Gaussian quadratures, and introduce methods for selecting sample points from a grid of phase-space localized Wilson basis. Both full configuration interaction and Hartree-Fock implementations are discussed and implemented numerically. As with many grid based methods, describing both tightly bound and diffuse orbitals is challenging so we have considered augmenting the Wilson basis set as projected Slater-type orbitals. We have also compared the Wilson basis set against the recently introduced wavelet transformed Gaussians (gausslets). Throughout, we give comments on the implementation and use small atoms and molecules to illustrate convergence properties of the Wilson basis.

Equations

$$w_{n}\left(x\right)=\left\{\begin{array}[]{ll}\phi\left(x-\frac{m}{2}\right)&\left[m\right]\in 2\,\mathbb{Z},k=0\\ \sqrt{2}\phi\left(x-\frac{m}{2}\right)\cos\left[2\pi kx\right]&\left[m\right]+k\in 2\,\mathbb{Z},k>0\\ \sqrt{2}\phi\left(x-\frac{m}{2}\right)\sin\left[2\pi kx\right]&\left[m\right]+k\in 2\,\mathbb{Z}+1,k>0\end{array}\right.,$$

$$\varphi \left(x\right)=\sum _{j,l\in \mathbb{Z}}{a}_{jl}\exp \left[2il\pi x\right]{\left(2v\right)}^{1/4}\exp \left[-v\pi (2x-j{)}^2\right].$$

$$d_{n}\left(x\right)=\left\{\begin{array}[]{ll}\frac{\eta_{n}}{\sqrt{{\scriptstyle\Delta}x}}\exp\left[{-\frac{\pi}{2{\scriptstyle\Delta}x^{2}}\left(x-m\,{\scriptstyle\Delta}x\right)^{2}}\right]&\left[m\right]\in 2\,\mathbb{Z},k=0\\ \frac{\eta_{n}}{\sqrt{{\scriptstyle\Delta}x}}\exp\left[{-\frac{\pi}{2{\scriptstyle\Delta}x^{2}}\left(x-m\,{\scriptstyle\Delta}x\right)^{2}}\right]\cos\left[\frac{\pi}{{\scriptstyle\Delta}x}k\left(x-m\,{\scriptstyle\Delta}x\right)\right]&\left[m\right]+k\in 2\,\mathbb{Z},k>0\\ \frac{\eta_{n}}{\sqrt{{\scriptstyle\Delta}x}}\exp\left[{-\frac{\pi}{2{\scriptstyle\Delta}x^{2}}\left(x-m\,{\scriptstyle\Delta}x\right)^{2}}\right]\sin\left[\frac{\pi}{{\scriptstyle\Delta}x}k\left(x-m\,{\scriptstyle\Delta}x\right)\right]&\left[m\right]+k\in 2\,\mathbb{Z}+1,k>0\end{array}\right.,$$

$${\eta }_n=2{\left(\frac{1}{2(1+(-1{)}^{\left[m\right]+k}\exp \left[-\pi k^2\right])}\right)}^{1/2}{\left(-i\right)}^{\left[m\right]+k}.$$

$$d_n\left(x\right)=$$

$$\left(\exp \left[i\frac{\pi }{\Delta x}k\left(x-m\Delta x\right)\right]+(-1{)}^{[m]+k}\exp \left[-i\frac{\pi }{\Delta x}k\left(x-m\Delta x\right)\right]\right)$$

$${\alpha }_n^{\pm }(x)=\frac{{\eta }_n}{\sqrt{\Delta x}}\exp \left[-\frac{\pi }{2\Delta x^2}{\left(x-m\Delta x\right)}^2\right]\exp \left[\pm k\frac{i\pi }{\Delta x}\left(x-m\Delta x\right)\right]$$

$$b_{n_1,n_2,...,n_D}(x_1,x_2,...,x_D)=\prod _{i=1}^D{d}_{n_i}\left(x_i\right)$$

$$S_{{\vec{n}}^{\prime }\vec{n}}=\prod _i{S}_{n_i^{\prime }{n}_i}=\int \mathrm{d}x{d}_{n_x}(x)d_{n_x^{\prime }}(x)\int \mathrm{d}y{d}_{n_y}(y)d_{n_y^{\prime }}(y)\int \mathrm{d}z{d}_{n_z}(z)d_{n_z^{\prime }}(z)$$

$$T_{{\vec{n}}^{\prime }\vec{n}}^{(j)}=T_{n_j^{\prime },n_j}\prod _{i\ne j}{S}_{n_i^{\prime }{n}_i}.$$

$$g\left(A,B,C,D\right)={\int }_{-\infty }^{\infty }\mathrm{d}xG\left(A,B,C,D,\mu ,x\right)$$

$$G\left(A,B,C,D,\mu ,x\right)=\exp \left[-A{\left(x-\mu \right)}^2\right]\exp [iB\left(x-\mu \right)]\exp [C]\exp [iD]$$

$$g\left(A,B,C,D\right)={\left(\frac{\pi }{A}\right)}^{1/2}\exp [-\frac{B^2}{4A}]\exp [C]\exp [iD]$$

$${\int }_{-\infty }^{\infty }\mathrm{d}xxG\left(A,B,C,D,\mu ,x\right)=\left(i\frac{B}{2A}+\mu \right)g(A,B,C,D)$$

$${\int }_{-\infty }^{\infty }\mathrm{d}x{x}^{2}G\left(A,B,C,D,\mu ,x\right)=\left(\frac{2A-B^2}{4A^2}+i\mu \frac{B}{A}+{\mu }^2\right)g(A,B,C,D)$$

$${\tilde{S}}_{n^{\prime }n}^{(a^{\prime }a)}$$

$$A$$

$$\begin{array}{rl}{\tilde{S}}_{n^{\prime }n}^{(a^{\prime }a)}& =\left({\eta }_n{\eta }_{n^{\prime }}\frac{1}{\Delta x}\right)g(A,B,C,D,\mu )\\ & ={\eta }_n{\eta }_{n^{\prime }}\exp \left[-\frac{\pi }{4}\left(m_{-}^2+{\left(k_{+}^{(a^{\prime }a)}\right)}^2\right)\right]\exp \left[-i\frac{\pi }{2}\left(m_{-}{k}_{-}^{(a^{\prime }a)}\right)\right].\end{array}$$

$$\begin{array}[]{ccl}S_{n^{\prime}n}^{(++)}&=&\tilde{S}_{n^{\prime}n}^{(++)}\\ S_{n^{\prime}n}^{(+-)}&=&(-1)^{\left[m\right]+k}\tilde{S}_{n^{\prime}n}^{(+-)}\\ S_{n^{\prime}n}^{(-+)}&=&(-1)^{\left[m^{\prime}\right]+k^{\prime}}\tilde{S}_{n^{\prime}n}^{(-+)}\\ S_{n^{\prime}n}^{(--)}&=&(-1)^{\left[m\right]+k+\left[m^{\prime}\right]+k^{\prime}}\tilde{S}_{n^{\prime}n}^{(--)}\end{array}$$

$$S_{n^{\prime }n}=\sum _a^{+,-}\sum _{a^{\prime }}^{+,-}{S}_{n^{\prime }n}^{(a^{\prime }a)}.$$

$$I_{n^{\prime }n}(O)=\sum _a^{+,-}\sum _{a^{\prime }}^{+,-}{I}_{n^{\prime }n}^{(a^{\prime }a)}(O)$$

$$I_{n^{\prime }n}^{(a^{\prime }a)}(x)$$

$$I_{n^{\prime }n}^{(a^{\prime }a)}(x^2)$$

$$\left(\frac{\mathrm{d}}{\mathrm{d}x}{\alpha }_{n^{\prime }}^{a^{\prime }}\right)\times \left(\frac{\mathrm{d}}{\mathrm{d}x}{\alpha }_n^a\right)$$

$$I_{n^{\prime }n}^{(a^{\prime }a)}\left({\frac{\mathrm{d}}{\mathrm{d}x}}^{†}\frac{\mathrm{d}}{\mathrm{d}x}\right)=\frac{{\pi }^2}{4\Delta x^2}\left(\frac{2}{\pi }-{\left(m_{-}+i{k}_{-}^{(a^{\prime }a)}\right)}^2\right)S_{n^{\prime }n}^{(a^{\prime }a)}.$$

$$T_{n^{\prime }n}=\sum _{a^{\prime },a}^{+,-}{I}_{n^{\prime }n}^{(a^{\prime }a)}\left({\frac{\mathrm{d}}{\mathrm{d}x}}^{†}\frac{\mathrm{d}}{\mathrm{d}x}\right)$$

$$\frac{1}{|r-r^{\prime }|}=\frac{2}{\sqrt{\pi }}{\int }_0^{\infty }\exp \left[-t^2(x-x^{\prime }{)}^2\right]\exp \left[-t^2(y-y^{\prime }{)}^2\right]\exp \left[-t^2(z-z^{\prime }{)}^2\right]\mathrm{d}t$$

$$I_{{\vec{n}}^{\prime }\vec{n}}^{{\vec{a}}^{\prime },\vec{a}}\left(\frac{1}{|r-r^{\prime }|}\right)={\int }_{-\infty }^{\infty }\frac{G_x(x)G_y(y)G_z(z)}{|r-r^{\prime }|}dxdydz$$

$${\tilde{V}}_{n_j^{\prime }{n}_j}^{(a_j^{\prime }{a}_j)}(t)$$

$$A_{jt}$$

$$C_{jt}$$

Full text

1. I Introduction

2. II Wilson Basis Functions

3. III Operator representation integrals

4. III.1 Overlap Integral

5. III.2 1D Kinetic energy operator

6. III.3 Coulomb integral

7. III.3.1 Electron-nuclear Coulomb integral

8. III.3.2 Two-electron Coulomb integral

9. IV Evaluating Coulomb matrix elements using Gaussian quadrature

10. V Calculating Energies

11. V.1 Choosing basis functions

12. VI One-Electron calculations

13. VI.1 Convergence of quadrature approximation

14. VI.2 Hermitian vs Product Approximation

15. VI.3 Using different sizes of modulated Gaussians phase-space localized functions

16. VI.4 Convergence 1S, 2S, and 2P hydrogen energy levels using different grid spacing

17. VII Two-Electron calculations

18. VIII Hartree-Fock

19. VIII.1 Four electron system using Hartree-Fock

20. IX Extensions of the Wilson basis set

21. IX.1 Using gausslets

22. IX.2 Augmenting the Wilson basis with Gaussian Slater type orbitals

23. IX.3 Other extensions

24. X Conclusion

Notation

${\scriptstyle\Delta}x$, Spacing between basis functions in position space. 2. 2.

$m$, Index for location of phase-space basis function in position space at $m{\scriptstyle\Delta}x$, not necessarily an integer but $m$ has integer spacing. 3. 3.

$k$, Index for location of phase-space basis function in momentum space. 4. 4.

$n$, Signifies phase-space basis function composite index $n={m,k}$. 5. 5.

$\vec{n}=\left[n_{1},n_{2},...n_{D}\right]$, The vector of all 1D indices for the multidimensional basis function with $D=3$ or $D=6$ 6. 6.

$L_{k}$, Number of grid functions in position space. 7. 7.

$L_{m}$, Number of grid functions in momentum space. 8. 8.

$L$, Total number of one dimensional phase-space functions. 9. 9.

$N$, Size of basis. 10. 10.

${{{\bf H}}}$, Matrices will be capital letters in bold and non-italicized. 11. 11.

${{{\bf v}}}$, Vectors will be small letters in bold and non-italicized. 12. 12.

$I(O)$, The integral for the operator $O$. 13. 13.

$G\left(A,B,C,D,\mu,x\right)$, The combination of two Gaussians. 14. 14.

$g\left(A,B,C,D\right)$, The integral of $G\left(A,B,C,D,\mu,x\right)$ over all $x$.

I Introduction

Computation of the electronic structure of a fixed nuclear potential occupies a swath of academic and commercial research in both the quantum and classical computational domains. In all cases, the basis set used to represent the possible locations of electrons plays a large role in the reach of a finite computational device. Both plane wave computations and spatially local basis functions form the major approaches. A key alternative, which includes the present study, makes use of basis functions localized in both momentum and real space. Such phase-space localized (PSL) basis functions in the form of wavelets have been used in electronic structure calculations with some success for a number of years.Arias (1999); Goedecker and Ivanov (1998); Harrison et al. (2004); Flad et al. (2002); Cho et al. (1993); Flad et al. (2006); Genovese et al. (2008); Shimshovitz and Tannor (2012a); White (2017) However, these basis functions are, with one exception,White (2017) defined only on a grid which makes obtaining analytic matrix elements very difficult for the unbounded Coulomb potential. In 1987, K. G. Wilson introduced a basis which is analytic everywhere and localized in phase-space.Wilson (1987); Sullivan et al. (1987)

LABEL:Sullivan1987 was the first to consider Wilson basis functions for electronic structure and suggests that PSL functions should have certain advantages over other grid basis methods. This early work only considered only a one-dimensional electronic problem and hints at the methodology for performing numerical calculations in real space. Here we develop the full machinery necessary to utilize Wilson basis sets very similar to the I. Daubechies, S. Jaffard, and J. JournDaubechies, Jaffard, and Journé (1991) construction for electronic structure calculations. The sum-of-products form of the full Hamiltonian in the Wilson basis allows the Hamiltonian to be decomposed into the separate Cartesian directions. This reduces memory requirements and can be used to perform matrix-vector products more efficiently using sequential summation.Wodraszka and Carrington (2017) We also use Gaussian quadrature to numerically evaluate the Hamiltonian terms. Throughout we use a variety of small electronic systems using both exact diagonalization (i.e. full configuration interaction) and the Hartree-Fock approximation. Further, our paper finds and describes reasonable numerical parameters and procedures for the Wilson basis and its extensions.

The organization of the paper is as follows. The Wilson basis is introduced in Section II. The representation of the operators in the Wilson basis is computed in Section III as closed form expressions. We then describe the evaluation of the Coulomb integrals using Gaussian quadrature in Section IV. Next, in Section V, we combine the integrals over the Wilson basis with sampling techniques to yield a complete computational procedure for calculating electronic energies. Then the parameters of algorithm are varied and tested on the exactly solvable one-electron hydrogen atomic system in Section VI. Section VII examines the Wilson basis applied to multiple states of the two-electron helium and molecular hydrogen systems. Section VIII discusses how Hartree-Fock calculations with the Wilson basis can be performed by using the methodology developed in the paper. In the penultimate section, Section IX, gausslet basis sets and standard Slater-type orbital basis sets are considered as possibilities to extend the Wilson basis. We summarize our conclusions and give an outlook in Section X.

Atomic units ($\hbar=m_{e}=a_{0}=4\pi\epsilon_{0}=1$) are used throughout.

II Wilson Basis Functions

The Wilson basis functions as defined in LABEL:Daubechies1991 are a product of 1D PSL functions that are of the form

[TABLE]

where modulated Gaussians are used as functions to generate $\phi\left(x\right)$ such that,

[TABLE]

In LABEL:Daubechies1991, the $a_{jl}$ are found either using a Zak transform or a convergent series in momentum space.

In order to simplify the calculations, we use the technique related to that of Refs. Poirier and Salam, 2004; Halverson and Poirier, 2012 but most similar to LABEL:Brown2016 to generate a Wilson basis. This involves using a grid of modulated Gaussians defined here as,

[TABLE]

where

[TABLE]

with ${\scriptstyle\Delta}x$ being the distance between functions in position space and $\left[m\right]$ is either the integer part of $m$ or the integer part of $m-1$, depending on where the $k=0$ are chosen to be placed. For these functions, we have $\Delta x\Delta p=\pi$. Each basis function is localized in position space around $m{\scriptstyle\Delta}x$ where the $L_{m}$ possible $m$ values have integer spacing but are in general real. The functions are also localized in momentum space at $\pm k\pi/{\scriptstyle\Delta}x$ where $k$ is a positive integer such that $k=0,1,2,3,...,L_{k}-1$. The index $n$ is taken to represent the composite index ${k,\left[m\right]}$ and the total 2D phase-space grid that represents one real-space dimension is composed of $L=L_{m}\times L_{k}$ functions. Each box in positive momentum space takes up $\pi$ but is combined with the corresponding negative momentum space partner such that the full basis function follows the uncertainty principle and is localized in $2\pi$ of phase-space.

As the modulated Gaussians of Eq. (3) are not orthogonal to each other, they must be orthogonalized in order to be used for calculating eigenvalues iteratively. We consider two schemes: symmetric orthogonalized ($\tilde{w}_{n}\left(x\right)=\sum_{j}{{{\bf S}}}_{nj}^{-1/2}d_{j}\left(x\right)$) and biorthogonalized ($b_{n}\left(x\right)=\sum_{j}{{{\bf S}}}_{nj}^{-1}d_{j}\left(x\right)$ combined with $d_{n^{\prime}}(x)$ as the dual). There are benefits to either choice but all resulting basis functions are localized exponentially in phase-space. The symmetric orthogonalized version ($\tilde{w}_{n}\left(x\right)$) of these basis functions are shown in Fig. 1 and form an orthonormal basis. When the basis functions $\tilde{w}_{n}\left(x\right)$ are not close to the boundary of the grid, the functions are very similar ($\max_{x}|\tilde{w}_{n}(x)-w_{n}(x)|\lesssim 10^{-3}$) to the I. Daubechies, S. Jaffard, and J. Journ (DJJ) functions $w_{n}\left(x\right)$. This can be seen in Fig. 2 for the $n=\left(0,1\right)$ and $n=\left(1,1\right)$ positions in phase-space. The two depicted functions are representative of all $m,k$ values in the interior of the phase-space domain. Equivalency was claimed in LABEL:Poirier2004a assuming the underlying modulated Gaussians used to generate the basis functions are the same. However, LABEL:Poirier2004a did not explicit show the equivalency and also did not appear to rigorously prove it. We do not make an effort to prove the equivalency of $w\left(x\right)$ and $\tilde{w}_{n}\left(x\right)$ but note that $\tilde{w}_{n}\left(x\right)$, $b_{n}\left(x\right)$ and $d_{n}\left(x\right)$ all have the properties of the original Wilson basis.

All functions ($\tilde{w}_{n}\left(x\right)$, $b_{n}\left(x\right)$ and $d_{n}\left(x\right)$) have symmetry properties in a checker board like fashion where $\left[m\right]+k\in 2\,\mathbb{Z}+1$ are odd and $\left[m\right]+k\in 2\,\mathbb{Z}$ are even. $\tilde{w}_{n}\left(x\right)$ and $b_{n}\left(x\right)$ are exponentially localized in phase-space while $d_{n}\left(x\right)$ is Gaussian localized. If a smaller grid of $d_{n}\left(x\right)$ is used to generate $\tilde{w}_{n}\left(x\right)$, exponential localization is retained but the functions are not as similar (especially for functions at the edge of the grid) to the construction of LABEL:Daubechies1991. That being said, the accuracy and pruneability of the basis is not impacted greatly.Halverson and Poirier (2012); Brown and Carrington (2016)

The symmetric (or anti-symmetric) position in momentum space is more explicitly obvious when Eq. (3) is in the complex exponential form with $\cos\left[px\right]=1/2\left(\exp\left[{{i\mkern 1.0mu}px}\right]+\exp\left[{-{i\mkern 1.0mu}px}\right]\right)$, $\sin\left[px\right]=-{i\mkern 1.0mu}/2\left(\exp\left[{{i\mkern 1.0mu}px}\right]-\exp\left[{-{i\mkern 1.0mu}px}\right]\right)$. This allows the Wilson basis of Eq. (3) to be written more succinctly as

[TABLE]

In Eq. (II), the values of $\left[m\right]\in 2\mathbb{Z}+1,k=0$ not included in the basis. If one ignored the normalization pre-factor, all $\left[m\right]\in 2\mathbb{Z}+1,k=0$ evaluate to zero.

We can further simplify notation by considering each the positive and negative momentum term separately in Eq. (II) using

[TABLE]

which results in $d_{n}(x)=\alpha_{n}^{+}(x)+(-1)^{\left[m\right]+k}\alpha_{n}^{-}(x)$.

To form a multidimensional basis, a product of the $1D$ basis functions of Eq. (II) for each of the $D$ dimensions is used. The number of dimensions is either three or six. In this section we are considering one electron integrals so $D=3$ but when considering the two electron integrals $D=6$. In both cases, the basis function is written as

[TABLE]

where there is a composite index $n_{i}$ for each of the $i=1,2,...,D$ dimensions.

III Operator representation integrals

This section outlines all the integrals needed to define the Hamiltonian representation in the Wilson basis. The one and two-body Coulomb integrals can be calculated advantageously as a sum-of-products. In electronic structure, the sum-of-product form has been previously used to assist in the evaluation of exchange and Coulomb integrals,Khoromskij, Khoromskaia, and Flad (2011) as well as evaluating the Møller-Plesset perturbation second order correction.Hohenstein, Parrish, and Martínez (2012) Exploiting grid based basis sets and a sum-of-product decomposition of the Hamiltonian into its three Cartesian product, has been applied to electronic structure, using a multi-resolution disjoint Legendre polynomial basis,Harrison et al. (2004) and a tensor decomposed sinc function basis.Jerke and Poirier (2018) Like the Wilson basis, the sinc functions are analytic. However, the Wilson basis functions are also localized in phase-space which has advantages for representing the cusps of electronic wavefunctions.

Due to the sum-of-products form, the overlap matrix and the Cartesian kinetic energy operator can be defined using a product of $1D$ matrix elements. All operator $O$ integrals will be performed using the notation of $I^{(a^{\prime}a)}_{n^{\prime}n}\left(O\right)$ which refers to the pair $\alpha_{n^{\prime}}^{a^{\prime}},\alpha_{n}^{a}$ with $a^{\prime},a\in\{+,-\}$ except for the overlap integral denoted $S^{(a^{\prime}a)}_{n^{\prime}n}$. The integral then includes four pairs of momentum combinations which are $(+k^{\prime},+k)$, $(+k^{\prime},-k)$, $(-k^{\prime},+k)$, and $(-k^{\prime},-k)$. The $I^{(+,+)}_{n^{\prime}n}\left(O\right)$ and $I^{(+,-)}_{n^{\prime}n}\left(O\right)$ integrals are equivalent to the complex conjugates of the $I^{(-,-)}_{n^{\prime}n}\left(O\right)$ and $I^{(-,+)}_{n^{\prime}n}\left(O\right)$ integrals respectively. The two-electron terms have a similar structure but require summing over more momentum combinations.

The overlap matrix for two basis functions in three dimensions is $\vec{n}^{\prime}=\left[n_{x}^{\prime},n_{y}^{\prime},n_{z}^{\prime}\right]$ and $\vec{n}=\left[n_{x},n_{y},n_{z}\right]$ is

[TABLE]

while the portion of the Cartesian Laplacian with second derivative in dimension $j$ is written as

[TABLE]

Throughout we will be using Gaussian integrals heavily so let us establish notation. We will always represent Gaussian integrals in the form

[TABLE]

where

[TABLE]

with $A,B,C,D,\mu$ are real constants. The integration of Eq. (10) results in

[TABLE]

The integrals of $x$ and $x^{2}$ are also presented as they are used to efficiently calculate the $\nabla^{2}$ operator below.

[TABLE]

and

[TABLE]

III.1 Overlap Integral

The overlap integral that needs to be evaluated for all four $(a^{\prime},a)$ combinations is given as

[TABLE]

where either $(-1)^{\left[m\right]+k}$ or $(-1)^{[m^{\prime}]+k^{\prime}}$ are excluded from $\tilde{S}^{(a^{\prime}a)}_{n^{\prime}n}$ but are included in the definition of $S_{n^{\prime}n}^{(a^{\prime}a)}$.

For Eq. (15), the appropriate values of $A,B,C,D$ for the product of $\alpha_{n}^{a}(x)$ and $\alpha_{n^{\prime}}^{a^{\prime}}(x)$ are

[TABLE]

where $m_{\pm}=m^{\prime}\pm m$ and $k^{(a^{\prime}a)}_{\pm}=a^{\prime}k^{\prime}\pm ak$. Thus, the overlap integral for a given $a^{\prime},a$ combination evaluates to

[TABLE]

Each of the $a^{\prime},a$ combinations are added together to form the full overlap matrix element (i.e. $S_{n^{\prime}n}=\tilde{S}_{n^{\prime}n}^{(++)}+(-1)^{\left[m\right]+k+\left[m^{\prime}\right]+k^{\prime}}\tilde{S}_{n^{\prime}n}^{(--)}+(-1)^{\left[m\right]+k}\tilde{S}_{n^{\prime}n}^{(+-)}+(-1)^{\left[m^{\prime}\right]+k^{\prime}}\tilde{S}_{n^{\prime}n}^{(-+)}$). The final definition of the partial overlap values is

[TABLE]

where the appropriate $(-1)^{\left[m\right]+k}$ or $(-1)^{\left[m\right]^{\prime}+k^{\prime}}$ are now included. The full overlap matrix for basis functions $d_{n}(x)$ with $d_{n^{\prime}}(x)$ is

[TABLE]

III.2 1D Kinetic energy operator

All 1D operators examined here can be written as

[TABLE]

where $I^{(a^{\prime}a)}_{n^{\prime}n}(O)$ are the partial integrals for operator $O$.

When evaluating the integrals associated with the kinetic energy operator, the values of (16) will be the same since the kinetic energy operator applied to a complex Gaussian results in the same Gaussian multiplied by a second degree polynomial. Therefore we will also need to integrate against the $x$ and $x^{2}$ operators to evaluate the Laplacian.

Using (13) and (14), the necessary partial integrals are given in terms of the overlap integral as

[TABLE]

and

[TABLE]

where $S^{(a^{\prime}a)}_{n^{\prime}n}$ are defined in Eq. (18) with $a,a^{\prime}\in\left\{+,-\right\}$.

To evaluate the integral of the $-\nabla^{2}$ operator, we will use the Cartesian representation $(-\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}-\frac{\mathrm{d}^{2}}{\mathrm{d}y^{2}}-\frac{\mathrm{d}^{2}}{\mathrm{d}z^{2}})$ and perform 1D integrals for $x,y,z$. For greater symmetry in the integrals, we will act with the derivative operator to the left and to the right with form $\frac{\mathrm{d}}{\mathrm{d}x}^{\dagger}\frac{\mathrm{d}}{\mathrm{d}x}\equiv-\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}$. The action of $\frac{\mathrm{d}}{\mathrm{d}x}^{\dagger}\frac{\mathrm{d}}{\mathrm{d}x}$ in the integral simplifies to,

[TABLE]

Using the $S^{(a^{\prime}a)}_{n^{\prime}n},I^{(a^{\prime}a)}_{n^{\prime}n}(x)$ and $I^{(a^{\prime}a)}_{n^{\prime}n}(x^{2})$ derived above and simplifying results in

[TABLE]

One can then use Eq. (20) to obtain the complete integral by summing over the four possible values of $a^{\prime},a$, to obtain

[TABLE]

which can be stored in an $L\times L$ matrix ${{{\bf T}}}$ and is the same for $x,y,z$.

III.3 Coulomb integral

The Coulomb integral is more challenging and will be the focus of the remainder of this section. We evaluate the integrals in closed form using the error function, $\textrm{erf}(x)=2\pi^{-1/2}\int_{0}^{x}\exp(-t^{2})dt$, as well as discuss the use of numerical quadrature to evaluate the functions quickly. The numerical integration is, in the end, preferred since it allows the integration to be done component-wise. As we will see below, the input to the error function has all Cartesian components combined.

To exploit the fact that our multi-dimensional basis functions, Eq. (7), are products of each Cartesian dimension, the form of the Coulomb operator we use is

[TABLE]

where $r=[x,y,z]$ and $r^{\prime}=[x^{\prime},y^{\prime},z^{\prime}]$. We will use Eq. (26), to perform the electron-nuclear and electron-electron integrations in close form.

III.3.1 Electron-nuclear Coulomb integral

The electron-nuclear integrals are performed using $r^{\prime}=r_{A}=[R_{x},R_{y},R_{z}]$ as the fixed nuclear position. We will evaluate the integral using the form Eq. (26); first by integrating over the one-dimensional coordinates, then by integrating over the dummy variable $t$ of Eq. (26). Since all terms including those introduced by the integral form of the Coulomb operator are all Gaussian, there are many simplifications along the way.

The required integral is given by

[TABLE]

The constants for $G_{j}(r)=G(A_{j},B_{j},C_{j},D_{j},\mu_{j},r)$ are the same as those given in Eq. (16). In general, these constant will depend on which direction is being discussed, hence the subscript.

We can simplify notation here and in subsequent sections by performing the integral over the individual spatial coordinates before simplifying the Gaussian expression which then remains parameterized by $t$. Consider for $j=x,y,z$ the integration

[TABLE]

The constants for the final Gaussian integral are given by

[TABLE]

Anticipating the final integration, we perform a change of variables to further convert the integral over $t\in(0,\infty)$ to $v\in(-1,1)$ with

[TABLE]

Performing the change of variables in Eq. (28) and integrating over $q$ results in

[TABLE]

where $c_{n_{i}^{\prime}n_{i}}^{(a_{i}^{\prime}a_{i})}=\frac{\pi}{8}(m_{i+}+{i\mkern 1.0mu}k_{+}^{(a^{\prime}a)}-2\frac{R_{i}}{{\scriptstyle\Delta}x})^{2}$. Our final definition for $V$ removes the prefactor

[TABLE]

The full 3D integral for a given $\vec{n},\vec{a},\vec{n}^{\prime},\vec{a}^{\prime}$ is then reduced to

[TABLE]

The complete integral for a nucleus at $r^{\prime}$ for a basis function $\vec{n}^{\prime}=\left[n_{x}^{\prime},n_{y}^{\prime},n_{z}^{\prime}\right]$ with $\vec{n}=\left[n_{x},n_{y},n_{z}\right]$ is then

[TABLE]

where the sum is performed over all $64$ combinations of $a_{x},a_{y},a_{z},a_{x}^{\prime},a_{y}^{\prime},a_{z}^{\prime}$ which can each be plus or minus. Evidently, the 3D Coulomb integral can be written in the multidimensional form of Eq. (20) with an additional integration over variable $v$. Eq. (34) can be evaluated in closed form to,

[TABLE]

III.3.2 Two-electron Coulomb integral

The form of the Coulomb operator for the two-electron integral is similar to Eq. (33) but defining $\vec{n_{i}}=\left[n_{x_{i}},n_{y_{i}},n_{z_{i}}\right]$ and $\vec{a_{1}}=\left[a_{x_{i}},a_{y_{i}},a_{z_{i}}\right]$ is,

[TABLE]

where

[TABLE]

Here $G_{i}^{k}=G(A_{i},B_{i},C_{i},D_{i},\mu_{i},q_{k})$ is the product of Gaussian corresponding to different dimensions $i=x,y,z$ and a one-dimensional electron coordinate $k=1,2$ along that that direction. To use this function we need to perform the integrals over $x_{1},y_{1},z_{1},x_{2},y_{2},z_{2}$ and then perform the integral over $t$. The change of variables here is $v=\frac{\pi-2{\scriptstyle\Delta}x^{2}t^{2}}{\pi+2{\scriptstyle\Delta}x^{2}t^{2}}$ which results in

[TABLE]

where

[TABLE]

where $b_{i}^{(a_{1}^{\prime},a_{1},a_{2}^{\prime},a_{2})}=\frac{\pi}{16}\left(m_{i_{1}+}-m_{i_{2}+}+{i\mkern 1.0mu}(k_{i_{1}+}^{(a_{i_{1}}^{\prime}a_{i_{1}})}-k_{i_{2}+}^{(a_{i_{2}}^{\prime}a_{i_{2}})})\right)^{2}$ with $i=x,y,z$.

The full 6-dimensional integral is given as a sum of $4^{6}=4096$ terms in closed form as

[TABLE]

Note that $\Delta x$ only appears as a prefactor to the integral. Thus, the Coulomb integrals do not need to be re-evaluated when changing the spatial distance between adjacent Wilson basis functions in real space.

IV Evaluating Coulomb matrix elements using Gaussian quadrature

In the previous section, we have evaluated all the overlap, kinetic and Coulomb integrals in closed form using the error function. However, the error function is not separable since the $x,y,z$ components are combined in the argument of the error function. In order to use a much larger number of basis functions, it is important to a) not store the full matrix b) increase the efficiency of matrix-vector products.

We can accomplish this using Gaussian quadrature over each dimension separately to get the Hamiltonian into a sum-of-products form. The overlap integral and the kinetic energy integrals are already a sum-of-products of one-dimensional operators so we only need to consider the Coulomb integrals.

Let us consider integrals of the form Eq. (34) which is in the form $\int_{-1}^{1}f(x)(1-x)^{a}(1+x)^{b}$ with $a=-0.5$ and $b=0$. Integrals of this form can be evaluated using Gauss-Jacobi quadrature Trefethen (2013) where the integral is approximated using $\sum^{N_{p}}_{\alpha}w_{\alpha}f(p_{\alpha})$. Here $N_{p}$ is the number of quadrature points $p_{\alpha}$ and their corresponding weights are $w_{\alpha}$.

This allows us to accurately convert the integral Eq. (33) into a sum-of-products form with $N_{p}$ terms such that,

[TABLE]

where

[TABLE]

with $i=x,y,z$. ${{{\bf V^{(1)}}}}(p_{\alpha},R_{i})$ is an $L\times L$ matrix with each element being a sum of $2^{2}=4$ terms.

Since Gaussian quadrature calculates polynomials of order $2N_{p}+1$ exactly, the convergence of the integral is superb. Test calculations have found that as few as seven quadrature points are necessary to obtain chemical accuracy as shown in the test results below. The same transformation can be performed for the two-body Coulomb as

[TABLE]

where

[TABLE]

with $i=x,y,z$. ${{{\bf V^{(2)}}}}(p_{\alpha})$ an $L^{2}\times L^{2}$ matrix with each element being a sum of $4^{2}=16$ terms. It is not possible to separate $x_{1},x_{2}$ in any obvious fashion analytically. This means that storing the two-electron matrix elements grows rapidly with the size of 1D sampling grid, and is responsible for most of the memory usage. The ${{{\bf V^{(2)}}}}(p_{\alpha})$ are the same for $x,y,z$ so only one matrix for each of the $N_{p}$ quadrature points need to be stored in memory.

The main value of these transformations is that we can now make use of the well established technique of performing sums sequentially, which is of common usage in the calculation of vibrational energy levels. The advantage is only realized for a “normal” (i.e. not generalized) eigenvalue problem. Using the modulated Gaussian grid of Eq. (II) directly results in a generalized eigenproblem unless we orthogonalize or biorthogonalize the representation. For the one-electron case, the Hamiltonian can be written in a sum-of-products form as

[TABLE]

where $O^{(t)}_{n_{i}^{\prime}n_{i}}$ is a 1D operator which is either the kinetic energy integral (of Eq. (25)), the overlap integral (of Eq. (19)) or one of the Coulomb terms (of Eq. (42)). The number of terms in the Hamiltonian representation is $N_{t}=3+N_{n}N_{p}$ with factor of three for the kinetic energy operator directions, $N_{n}$ is the number of nuclei, and $N_{p}$ is the number of quadrature points used to evaluate each nuclei-electron Coulomb term. Each ${{{\bf O_{i}^{(t)}}}}$ for dimension $i$ is of size $L\times L$. This means that the full generalized eigenvalue problem can be written as

[TABLE]

where ${{{\bf Z}}}$ is the matrix of eigenvectors, and ${{{\bf E}}}$ is the diagonal matrix of eigenvalues. The total memory requirements to store the sum-of-products matrix representation is $3N_{p}L^{2}$ which is much less than the full matrix which requires $L^{6}$. The ${{{\bf Z}}}$ and ${{{\bf E}}}$ are of size $L^{3}\times L^{3}$. Eq. (46) can easily be converted to an eigenvalue problem with no overlap matrix by either using

[TABLE]

with ${{{\bf U}}}=({{{\bf S}}}^{1/2}\otimes{{{\bf S}}}^{1/2}\otimes{{{\bf S}}}^{1/2}){{{\bf Z}}}$ or

[TABLE]

with ${{{\bf U}}}=({{{\bf S}}}\otimes{{{\bf S}}}\otimes{{{\bf S}}}){{{\bf Z}}}$. Both forms have been used successfully in previous vibrational calculations.Brown and Carrington (2015a) Eq. (47) has the advantage of being Hermitian while Eq. (48) generally produces a smaller basis representation. Most calculations in this study use Eq. (47) but both are examined. We also emphasize that Eq. (47) will result in a calculation very similar to using the Wilson basis of LABEL:Daubechies1991 if a large grid of Eq. (II) is used to generate the ${{{\bf O_{i}}}}^{(t)}$ and ${{{\bf S}}}$ matrices.

Similarly, the two-body Coulomb term can be orthogonalized as

[TABLE]

or biorthogonalized as

[TABLE]

where ${{{\bf O}}}^{(t)}_{x_{1},x_{2}}$ represents the full 4D matrix with $V^{(2)}_{n_{i_{1}}^{\prime}n_{i_{1}}n_{i_{2}}^{\prime}n_{i_{2}}}$ as elements and ${{{\bf S}}}_{1,2}={{{\bf S}}}\otimes{{{\bf S}}}$. It is clear that $i_{1},i_{2}$ are not separable but we can do each pair $\{x_{1},x_{2}\},\{y_{1},y_{2}\},\{z_{1},z_{2}\}$ separately.

V Calculating Energies

In order to calculate eigenvalues efficiently, it is important to perform matrix-vector products sequentially. In order to simplify the notation, we restrict this attention to the one-electron problem and define ${{{\bf B}}}_{j}^{(t)}={{{\bf S}}}^{-1/2}{{{\bf O}}}^{(t)}_{x}{{{\bf S}}}^{-1/2}$ or ${{{\bf O}}}^{(t)}_{x}{{{\bf S}}}^{-1}$ such that the Hamiltonian is now

[TABLE]

where the first summation is over all the terms and $t$ is an arbitrary labelling of the term. The full matrix vector product for one coefficient $\vec{u}_{\vec{n}^{\prime}}$ corresponding to indices $\vec{n}^{\prime}=[n_{x}^{\prime},n_{y}^{\prime},n_{z}^{\prime}]$ is given as

[TABLE]

with the sum performed over all $L$ values on the grid for each of $x,y,z$. The total summation performed naively using Eq. (52) requires a scaling of $L^{6}$ since, for each of the $L^{3}$ components of ${{{\bf w}}}$, one must sums over all $L^{3}$ components of ${{{\bf u}}}$. One can instead perform matrix-vector products sequentially using Eq. (53). This is performed by using intermediate vectors $\vec{{{{\bf u}}}}^{\prime}$ and $\vec{{{{\bf u}}}}^{\prime\prime}$. First, $\vec{{{{\bf u}}}}^{\prime}$ is formed by performing the sum over $z$ first such that,

[TABLE]

which involves a sum over $L$ terms for $L^{3}$ indices. Then $\vec{{{{\bf u}}}}^{\prime\prime}$ is generated by taking $\vec{{{{\bf u}}}}^{\prime}$ as an input for the sum over $y$ such that,

[TABLE]

which involves another sum over $L$ terms for each $L^{3}$ indices. The sum is then performed over $x$ using $\vec{{{{\bf u}}}}^{\prime\prime}$ as an input which results in

[TABLE]

Two advantages can be noted. First, the three step process of Eq. (54)-Eq. (56) requires $3L^{4}$ operations compared to $L^{6}$ for Eq. (52) with the resulting vectors $\vec{u}_{\vec{n}^{\prime}}$ being equivalent. This is the computational benefit of having a sum-of-products Hamiltonian. The other key benefit is memory usage. One only needs to store a number of $L\times L$ matrices compared to the full Hamiltonian of $L^{3}\times L^{3}$ size.

For the two-electron problem, the matrices are not separable for $j_{1},j_{2}\in\{x,y,z\}$. This means that the sum that needs to be performed is

[TABLE]

where $V_{n_{j_{1}}^{\prime}n_{j_{1}}n_{j_{2}}^{\prime}n_{j_{2}}}$ are the $L^{2}\times L^{2}$ matrices for the two-electron Coulomb operator. Performing this matrix vector product sequentially requires $3\times L^{6}\times L^{2}=3L^{8}$ operations. If the Coulomb term was not separated into $3$ products, the total cost of the matrix-vector product would be $L^{12}$. Storing these $L^{2}\times L^{2}$ matrices is the major RAM requirement of these calculations since a fairly large grid is required to converge eigenvalues. Storing the full matrix would not be feasible.

The direct product basis for a two-electron problem has $L^{6}$ components but even this becomes intractable quickly. Therefore, instead of using the full direct product grid, we only include basis functions that have significant overlap with the desired wavefunction(s). We call this set of basis function labels, $\vec{n}$, as $\beta$. The size of $\beta$ is the number of one-electron basis functions denoted $M$. Using a general pruned basis $\beta$ complicates the matrix-vector products significantly. The matrix representation is now

[TABLE]

where ${{{\bf P}}}_{\beta}$ is a rectangular matrix that projects out the appropriate matrix elements for basis $\beta$. The sequential summation is now,

[TABLE]

where $n_{z}(n_{x},n_{y})$ denotes that the summation for dimension $z$ is performed only over those values of $n_{z}$ that are in the set $\beta$ that have corresponding $n_{y},n_{z}$ values. This is the effect of the projector ${{{\bf P}}}_{\beta}$. The second summation is performed over all values of $n_{y}(n_{x})$ which indicates that only $n_{y}$ values that have a corresponding $n_{x}$ are included while $n_{z}$ can be any value in $\beta$. The final summation is performed over all $n_{x}$ in $\beta$ but only generate the output for $\vec{n}^{\prime}\in\beta$ which is the effect of the projector ${{{\bf P}}}^{T}_{\beta}$. This means that in order to perform matrix-vector products sequentially, intermediate vectors $\vec{{{{\bf u}}}}^{\prime},\vec{{{{\bf u}}}}^{\prime\prime}$ will have a larger size than $M$.

If one uses the form of Eq. (59) as stated, then the intermediate vectors $\vec{{{{\bf u}}}}^{\prime},\vec{{{{\bf u}}}}^{\prime\prime}$ will most likely be larger than necessary. One only needs to retain any intermediate basis functions that includes the union of set $\left\{n_{x},n_{y},\bar{n_{z}}\right\}$ and $\left\{\bar{n_{x}},n_{y},n_{z}\right\}$, with the $\bar{n_{i}}$ signifying that all values of $n_{i}$ are taken while only the combinations $n_{j},n_{k}$ in the retained basis $\beta$ are included. Therefore, the intermediate vectors will still be greater than $M$ but smaller than what would result using Eq. (59). This is explained well in LABEL:Wodraszka2017 for any number of dimensions and is implemented here.

There is also an approximation that can be made such that intermediate vectors are always the same size. This is the product approximation where the matrix-representation of Eq. (47) is replaced by,

[TABLE]

where ${{{\bf F}}}_{x}={{{\bf P}}}_{\beta}^{T}\left({{{\bf S}}}^{-1/2}{{{\bf O}}}^{(t)}_{x}{{{\bf S}}}^{-1/2}\otimes{{{\bf I}}}_{y}\otimes{{{\bf I}}}_{z}\right){{{\bf P}}}_{\beta}$, ${{{\bf F}}}_{y}={{{\bf P}}}_{\beta}^{T}\left({{{\bf I}}}_{x}\otimes{{{\bf S}}}^{-1/2}{{{\bf O}}}^{(t)}_{y}{{{\bf S}}}^{-1/2}\otimes{{{\bf I}}}_{z}\right){{{\bf P}}}_{\beta}$, and ${{{\bf F}}}_{z}={{{\bf P}}}_{\beta}^{T}\left({{{\bf I}}}_{x}\otimes{{{\bf I}}}_{y}\otimes{{{\bf S}}}^{-1/2}{{{\bf O}}}^{(t)}_{z}{{{\bf S}}}^{-1/2}\right){{{\bf P}}}_{\beta}$ where ${{{\bf I}}}_{i}$ is the identity operator for coordinate $i$. Performing matrix-vector products is now

[TABLE]

The main disadvantage of the product approximation is that the Hamiltonian representation is no longer Hermitian. To justify this claim, consider an approximate sum-of-products approximation where $H_{trunc.}=\sum_{i}h_{i}h_{i+1}$ with Hermitian conjugate $H_{trunc.}^{\dagger}\approx\sum_{j}h_{i+1}h_{i}$. The action of the operator and its dual are only the same when $h_{i}h_{i+1}=h_{i+1}h_{i}$. This can be rectified by taking the transpose of the original ordering of the above matrices such that ${{{\bf F}}}_{x}{{{\bf F}}}_{y}{{{\bf F}}}_{z}$ becomes $\frac{1}{2}({{{\bf F}}}_{x}{{{\bf F}}}_{y}{{{\bf F}}}_{z}+{{{\bf F}}}_{z}{{{\bf F}}}_{y}{{{\bf F}}}_{x})$. See Ref. Cooper and Carrington (2009) for the first numerical application to molecular physics. However, in our context, this approach doubles the computational cost and as will be shown later in Section VI, is not necessary for accurate calculations to be made. One simply needs to use Arnoldi iterations as opposed to Lanczos iterations to calculate eigenvalues/eigenvectors Lehoucq, Sorensen, and Yang (1998). This does increase memory requirements as a set of basis vectors need to be stored but, for the calculations performed here, the main memory costs are from storing the Hamiltonian matrix elements for the two-electron Coulomb terms.

V.1 Choosing basis functions

It is impossible to know a priori the overlap of the basis functions $d_{n}(x)$ with the desired eigenfunctions. Choosing functions that are centered in the classically allowed region can provide semi-quantitative accuracy,Halverson and Poirier (2012) however this is difficult to do in multiple dimensions. A better method is to iteratively improve the basis function by starting with a small basis (motivated by classical phase-space energies) and progressively add more functions around the most important basis functions to improve the description of the wavefunction.Pipek and Nagy (2012); Shimshovitz, Bacic, and Tannor (2014); Brown and Carrington (2015b) In this study, the functions included in set $\beta$ are determined by considering the diagonal elements of the density matrix for the eigenstate.

Consider first the case where we are interested in the ground state and we will generalize to cases where we consider more than one state subsequently. At the start of the computation, $L$ is fixed and $\beta$ is a subset of $M$ functions chosen from the full $L^{3}$ basis functions. After each computation of an approximate wave function, we have $v_{\vec{n}i}$ as the coefficient for basis function $\alpha_{\vec{n}}$. The importance function of each basis function is then defined by $\mathcal{P}_{\vec{n}}=v_{\vec{n}i}^{2}$. If one is interested in optimizing the pruned basis for more than one eigenstate at a time, then the total importance function is the sum of the importance vector for each targeted state.

The importance vector determines which basis function are pruned and determines which basis functions will be added for the next iteration. Those basis function with importance values below a cutoff threshold are removed and not considered in future iterations. Expansion of the set $\beta$ occurs near the functions with the largest importance values. Additions are made in all $3\times 4$ directions in phase-space: three for each dimension $x,y,z$ and four for each choice of $\pm m,\pm k$. Special care needs to be (and was) taken to exclude redundant function that are already in $\beta$ and to exclude $\left[m\right]\in 2\mathbb{Z}+1,k=0$. For basis functions of $m=l,k=0$ the additions to the set $\beta$ are only $m=l\pm 2,k=0$, $m=l\mp 1,k=1$ and $m=l,k=1$.

For all calculations in this paper, we are using basis functions in $x,y,z$ that can be symmetry adapted with respect to inversion of single coordinate e.g. $f(-x,y,z)=\pm f(x,y,z)$. All systems studied in this paper have symmetry about the origin in at least two coordinates and therefore fewer basis functions will be needed to represent their wavefunction. This is true even though the full wavefunction (including correlation) does not respect this single-coordinate inversion symmetry. The asymmetry introduced by correlation is a much smaller portion of the total wavefunction and therefore requires fewer anti-symmetric basis functions to describe it accurately. This means that at the origin, $k\in 2\mathbb{Z}$ are symmetric functions while $k\in 2\mathbb{Z}+1$ are anti-symmetric functions. All other functions are localized in four positions in phase-space ($\pm m{\scriptstyle\Delta}x,\pm k\pi/{\scriptstyle\Delta}x$) except for $k=0$ which is localized at two points. The indices are then relabeled such that positive $m$ means moving outward from the origin in symmetric functions while a decrease in negative $m$ from means moving outward from the origin in anti-symmetric functions.

VI One-Electron calculations

There are four parameters that one controls, the lattice spacing in position space ${\scriptstyle\Delta}x$, the number (and type) of quadrature points used for the Coulomb terms, the number of underlying $d_{n}\left(x\right)$ functions, and whether the Hermitian (Eq. (59)) or product approximated (Eq. (61)) Hamiltonian representation is used. The hydrogen atom is used as the test for all parameters. The basis is symmetry adapted for $x,y,z$, with basis functions located at $...,-2{\scriptstyle\Delta}x,-{\scriptstyle\Delta}x,0,{\scriptstyle\Delta}x,2{\scriptstyle\Delta}x,...$ and the hydrogen nucleus located at $0,0,0$. The size of the position grid is always of the size $L_{k}=4\mathbb{Z}-1$ so that $m=0,k=0$ is in the basis with with an odd number of $m$ values on either side. The symmetric orthogonalized form of Eq. (47) with Hermitian matrix-vector products is used unless otherwise stated. Also, the underlying $d_{n}\left(x\right)$ grid is only has large as necessary except for section VI.3. Implicitly restarted Arnoldi via ARPACKLehoucq, Sorensen, and Yang (1998) is used to calculate eigenvalues/eigenvectors. When performing test of grid spacing and when comparing the Hermitian and product approximations, $19$ Gauss-Jacobi quadrature points were used. However, 13 quadrature points were used for all multi-electron calculations consistent with our findings in the next subsection.

VI.1 Convergence of quadrature approximation

The convergence of the Gauss-Jacobi quadrature approximation can be seen in Fig. 3 and Fig. 4 for the 1S (with ${\scriptstyle\Delta}x=\sqrt{\pi}$) and the 2P state (with ${\scriptstyle\Delta}x=\sqrt{\pi}$) respectively. For the 1S state with ${\scriptstyle\Delta}x=\sqrt{\pi}$, chemical accuracy $\epsilon<10^{-3}$ Hartree can be achieved with only $7$ quadrature points. This is not true if one uses ${\scriptstyle\Delta}x=2\sqrt{\pi}$ where $11$ quadrature points are required. The reason for this is that the maximum momentum index required for chemical accuracy is $k=15$ with ${\scriptstyle\Delta}x=2\sqrt{\pi}$ while only $k=7$ for ${\scriptstyle\Delta}x=\sqrt{\pi}$. A more highly oscillatory basis requires a higher degree polynomial to represent and therefore more quadrature points.

Therefore, all later calculations are performed with 13 Gauss-Jacobi quadrature points which is more than enough to achieve chemical accuracy without the quadrature error entering the calculation. One can also use Gauss-Legendre quadrature but the number of points required to achieve chemical accuracy is generally larger by about $50\%$.

VI.2 Hermitian vs Product Approximation

For the 1S state of hydrogen with a position spacing ${\scriptstyle\Delta}x=\sqrt{\pi}$, the comparison between the Hermitian (of Eq. (59)) and product approximated form (of Eq. (61)) is shown in Fig. 5. It is clear that the product approximation only manifests when high accuracy is required. For this calculation, that occurs after chemical accuracy has been achieved. This is not the case if ${\scriptstyle\Delta}x=2\sqrt{\pi}$ for which the difference is evident when the error is approximately $0.002$ $E_{h}$ and error below $0.001$ $E_{h}$ requires a larger basis. Even with this disadvantage, the speed of the calculation using the product representation is much faster as intermediate basis sizes can be upwards of six times larger than the retained basis.

VI.3 Using different sizes of modulated Gaussians phase-space localized functions

We have also testing using different sizes of modulated Gaussian (Eq. (II)) grids for the 1S state of a hydrogen atom. The large grid includes indices $m\in\left[-21,21\right],k\in\left[0,22\right]$. The small grid has indices $m\in\left[-7,7\right],k\in\left[0,13\right]$. The basis functions after the 40 expansion iterations shown have indices $m\leq 6$ and $k\leq 11$. In Fig. 6, the difference between the large and small grid of $d_{n}\left(x\right)$ functions are insignificant near convergence. In fact, the coefficients of the eigenvector are the same to several (approximately four) digits for a given $n$. Knowing this, we only use as many functions as is necessary to cover a large enough region of phase-space to converge the calculations. This is especially important for the two-electron calculations when memory resources become the constraining factor.

VI.4 Convergence 1S, 2S, and 2P hydrogen energy levels using different grid spacing

The S orbitals are determined with only even basis functions for $x,y,z$ while the 2P state is determined by using even basis functions for $x,y$ and odd basis functions for $z$. Without symmetry adaption, the number of basis functions would be $8=2^{3}$ times larger.

The convergence of hydrogen energy levels is determined by growing the basis using only the wavefunction from the desired state (see Subsection V.1). The convergence has been examined for the spacing of ${\scriptstyle\Delta}x=\sqrt{\pi}$ and ${\scriptstyle\Delta}x=2\sqrt{\pi}$ with each state being optimized independently. The discontinuous derivative of the 1S and 2S states results in a slower convergence than for the 2P state. For ${\scriptstyle\Delta}x=\sqrt{\pi}$ (shown in Fig. 7), the convergence of the 1S, 2S and 2P states is $N^{1.63}$, $N^{1.72}$ and $N^{2.236}$ respectively. When ${\scriptstyle\Delta}x=2\sqrt{\pi}$, the $2S$ and $2P$ states converge much more quickly ($N^{2.602}$ and $N^{5.563}$ respectively) compared to $1S$ only having a convergence of $N^{0.901}$. This is why multiresolution wavelet PSL basis sets are desirable White (2017). Unfortunately, the Wilson basis does not have this property but we return to this point in Section IX.1.

VII Two-Electron calculations

To perform two-electron calculations, we simply take the set of one electron basis functions and make a direct-product to add the other three dimensions.

[TABLE]

where $\vec{n}_{1}=[n_{1},n_{2},n_{3}]$ and $\vec{n}_{2}=[n_{4},n_{5},n_{6}]$ are both contained in the pruned basis function set $\beta$. This means that no product approximation is made in the individual two-electron Coulomb terms for $x,y,z$ regardless of whether one uses Eq. (59) or Eq. (61). The benefit of using Eq. (62) is that spin can be taken into account if desired.

For example, one can differentiate between singlet and triplet states by enforcing symmetry under exchange. This is explicitly written as,

[TABLE]

for the singlet state and

[TABLE]

for the triplet state. Even when using the product approximation (and Arnoldi iterations), symmetry adaptation isolates singlet and triplet states effectively as the product approximation error is small.

If one uses Eq. (59) then singlet/triplet labels can be obtained with a single matrix-vector product for each Lanczos iteration using symmetry adapted Lanczos.Wang and Carrington (2001) This is because the symmetry projection operator ${{{\bf P}}}$ commutes with the Hamiltonian such that ${{{\bf PHv}}}={{{\bf HPv}}}$. Therefore, a sum of vectors with different symmetries applied to the Hamiltonian, followed by the symmetries then being projected out is equivalent to applying the Hamiltonian to each vector with distinct symmetry separately. One then stores, and calculates eigenvalues from, separate Lanczos tridiagonal matrices for each symmetry. Note that symmetry adapted Lanczos was not implemented in this pilot study.

The first test case is for the ${}^{1}S$ ground state of the helium atom. The convergence of the basis with ${\scriptstyle\Delta}x=0.45\sqrt{\pi}$ on a grid of $L_{k}=11,L_{m}=11$ is found to be $M^{-0.726}$ (where $M=N^{2}$) for the symmetric basis and $N^{-0.781}$ for the biorthogonal basis as shown in Fig. 9. The most accurate energy calculated is found to be $-2.9013$ $E_{h}$ with the exact solution being $-2.90372...$ $E_{h}$Kurokawa, Nakashima, and Nakatsuji (2008). This required a basis of $3,325\times 3,325=11,055,625$ functions. This calculation took just under five hours using 16 cores on a Xeon(R) E5-2640 processor and required $28.5$GB of RAM. The majority of the RAM usage ($22$ GBs) were required to store the two-body Coulomb matrix representation at each of the 13 quadrature points. A value of $-2.9017$ $E_{h}$ using the biorthogonal representation used $2603\times 2603=6,775,609$ basis functions and required just under $3$ hours of computing time using $26.2$GB of RAM on the same computer. The reduction of approximately $20\%$ in the one-electron basis size is fairly consistent across all calculations examined when using the biorthogonal representation relative to the orthogonal basis.

We also performed calculations on the hydrogen molecule at a distance of $0.74$Å  and obtained the lowest two singlet states along with the lowest energy triplet state. The sampling grid was $L_{m}=11,L_{k}=11$ with ${\scriptstyle\Delta}x=1.3\sqrt{\pi}$. Optimizing the ground state separately resulted in an energy of $-1.8872$ $E_{h}$ with a basis size of $3180\times 3180=10,112,400$. When optimizing for all three eigenvalues, a basis of $2898\times 2898=8,398,404$ obtained energies of $-1.8865$ $E_{h}$ for $X\,^{1}\Sigma_{g}^{+}$, $-1.4961$ $E_{h}$ for $b\,^{3}\Sigma_{u}^{+}$ and $-1.4254$ $E_{h}$ for $B\,^{1}\Sigma_{u}^{+}$. All three of these values are more converged than those calculated with LABEL:Jerke2018. This is especially true for $B\,^{1}\Sigma_{u}^{+}$ which is $0.007E_{h}$ lower in energy.

VIII Hartree-Fock

Hartree-Fock or mean-field ansatz is also possible using the Wilson basis functions. This is because the Coulomb basis is not truncated but rather remains a direct product of basis functions, Eq. (62). The Hartree-Fock equations are derived by taking the gradient of the energy functional of charge density matrices. The resulting operator is the Fock matrix given in standard notation as Szabo and Ostlund (1996)

[TABLE]

with $(\vec{n}_{1}\vec{n}_{1}^{\prime}|\vec{n}_{2}\vec{n}_{2}^{\prime})=I_{\vec{n}_{1}^{\prime}\vec{n}_{2}^{\prime},\vec{n}_{1}\vec{n}_{2}}(|r_{1}-r_{2}|^{-1})$. We can achieve a more efficient algorithm by rearranging the action of the second term as follows

[TABLE]

This second form allows us to perform the action of the Fock operator without constructing the matrix. In some cases, the first form could be more efficent if the number of electrons is very high but the second form is what has been implemented. The integral driven procedure scales as $3N_{e}(N^{2}+N^{8/3})/2$. If instead one first forms the charge density matrix, $P_{\lambda\sigma}=\sum_{a}^{N_{e}/2}C_{\lambda a}C_{\sigma a}$, then performs the remaining summations in Eq. (66), a cost of $(N^{2})(3N^{4/3})/2$ is expected. It is therefore cheaper to use the density matrix approach only when the number of electrons is less than $N^{2/3}$. For all calculations performed here, $N_{e}\ll N^{2/3}$, so the integral driven approach was utilized.

The convergence of the Hartree-Fock calculation of the helium atom is shown in Fig. 10. The convergence using the symmetric basis is $N^{-1.561}$ which is essentially the same convergence as the full two electron calculation, but $N$ is denoted by the size of the one-electron basis here. The most accurate value obtained is $-2.86012$ $E_{h}$ using a basis of $2525$ functions. The exact value is $-2.8618$ for an error of less than $0.002$ $E_{h}$.

VIII.1 Four electron system using Hartree-Fock

Using Eq. (67), we study the four-electron system LiH at a spacing of $3$ $a_{0}$ with $H$ at $-0.75\,a_{0}$ and Li at 2.25$a_{0}$. The system geometry is such that both atoms are on the $x$-axis. Consequently, only even-parity functions are necessary for $y$ and $z$ directions reducing the necessary basis size by a factor of four. By choosing the step size in each direction the Couloumb matricies are the same in each direction. By performing the same transform to the symmetry adapted basis to each coordinates maintains the equivalence of the Coulomb matrices in each direction. Hence, we can store a single $L\times L\times L\times L$ Couloumb matrix at each quadrature point rather than maintaining separate representations in each Cartesian direction.

The basis functions for this calculation were sampled from a grid of $L_{m}=15$, $L_{k}=7$ with $N_{p}=13$ quadrature points and required $\approx 17$GB of memory. The majority of this memory overhead comes from the thirteen $(15\times 7)^{4}$ two-electron matrices representing approximately 13 GB of memory.

To test the performance of the Wilson basis on the LiH system, we compared against the standard basis sets of STO-3G and 3-21G with respective Hartree-Fock energies of $-7.8623$ and $-7.9295$ respectively. Of the total $L^{3}=1,157,625$ Wilson basis functions, a pruned basis set $\beta$ of size $N=2048$ was necessary to achieve STO-3G accuracy and $N=3511$ was needed to reach 3-21G accuracy.

We performed the calculations using 16 cores of a Xeon(R) E5-2640 processor with 64GB RAM, part of Dartmouth’s Discovery Cluster. To obtain STO-3G accuracy took 10.3 hours while obtaining 3-21G accuracy took approxmately $38$ hours.

IX Extensions of the Wilson basis set

The necessity of describing both diffuse and tightly bound eigenfunctions simultaneously complicates the usage of Wilson basis functions. Therefore, the Wilson basis will probably be more successful when combined with pseudopotentials. Instead of using pseudopotentials, two alternative ideas are tested.

The first is to utilize a multiresolution basis and the second is to augment the Wilson basis with projected Slater type orbitals. The multiresolution gaussletsWhite (2017) are more effective at describing both diffuse and tightly bound eigenfunctions but, as we point out in the next subsection, there are cases where the Wilson basis remains better.

The second method makes use of the sum-of-products form of the STO-$n$G basis. We showed that accurate energies of one electron in a LiH2 potential can be obtained by replacing high momentum functions with the STO-$6$G basis centered on the Li nuclei. This second technique will likely be most useful when using an expansion of the eigenfunction as a sum-of-products described in LABEL:Thomas2018 and LABEL:Jerke2018.

We will close this section with a tour through other ideas for improving the overall numerical implementation. The final subsection will highlight several paths for improvement that were not done in this study.

IX.1 Using gausslets

Earlier in section VI.4, we noted that the convergence rates of the 1S state is faster when the the Wilson function grid spacing is small but the 2S and 2P states converge faster with larger grid spacing. This suggests that a single choice of resolution will not work for multiple states. In this subsection, we explore multiresolution approaches and compare them against the Wilson basis. We find that the convergence rates of the Wilson basis are faster for individual states but the multiresolution functions can better capture simultaneous convergence of multiple states.

It has been shownShimshovitz and Tannor (2012a) that using multiresolution wavelets for Coulomb potentials can provide a significant reduction in the number of basis functions required for one-dimensional convergence. However, most wavelets used previously are only defined on a grid, and as the electron-nuclei and electron-electron potentials are unbounded, the Hamiltonian representations become non-trivial when performing quadrature due to the singulaity at various points in space. One then either needs to use a pseudo-potentialGenovese et al. (2008) or soften the discontinuityShimshovitz and Tannor (2012a).

In another approach to this problem, WhiteWhite (2017) developed a wavelet theory using “gausslets” that are defined as certain linear combinations of a equally spaced grid of Gaussians.

This allows a multi-resolution wavelet transform to be performed on a grid of Gaussian functions. Starting from a single wavelet, a set of three transformations can be used to decrease the resolution systematically: a transformation to an even function, a transformation to an odd, and a scaling transform. The even/odd symmetries are defined by inversion about the center of the gausslet. Here we consider a calculation for the electronic energy of the hydrogen atom.

The grid of Gaussians is defined as

[TABLE]

with $\delta x$ between each Gaussian on the grid. For each cartesisan direction, we consider a one-dimensional grid with $i$ ranging between $-4208$ and $4208$. The overlap integral for two Gaussians on the grid is

[TABLE]

The kinetic energy operator integral is

[TABLE]

and the one-electron Coulomb integral for each Legendre quadrature point $p_{\alpha}$ is

[TABLE]

The full Coulomb operator is

[TABLE]

We use the $\mathcal{G}_{10}$ gausslet with application of the $W_{652}$ wavelet transform repeated $4$ times and a Gaussian grid spacing of $\delta x=0.015$. This choice was found to give the best results using the methodology described in the present paper. The matrix elements were calculated by generating the transformation matrix ${{{\bf G}}}$ for each basis function at position $n=(m,k)$

[TABLE]

with $G_{ni}$ representing a vector of coefficients in ${{{\bf G}}}$ for $n=(m,k)$ such that $m=-8,-7..,0,...,7,8$ for $k=0$ and $m=-7.5,...-0.5,0.5,...,7.5$ for the even and odd wavelets at $k=1,2,...,8$. The centers of the wavelets change depending on the $k$-dependant spacing of the wavelets ${\scriptstyle\Delta}x^{(k)}$ given by ${\scriptstyle\Delta}x^{(7)}={\scriptstyle\Delta}x^{(8)}=9\delta x$, $\,{\scriptstyle\Delta}x^{(5)}={\scriptstyle\Delta}x^{(6)}=27\delta x$, $\,{\scriptstyle\Delta}x^{(3)}={\scriptstyle\Delta}x^{(4)}=81\delta x$ and $\,{\scriptstyle\Delta}x^{(0)}={\scriptstyle\Delta}x^{(1)}={\scriptstyle\Delta}x^{(6)}=243\delta x$. The even and odd wavelet transforms shift the center by a factor of ${\scriptstyle\Delta}x_{k}/2$ and is responsible for the different $m$ indexing between $k\neq 0$ and $k=0$.

We generate the basis by starting with gausslet $\mathcal{G}_{10}$,White (2017) defined by the coefficient vector $\vec{G}^{(10)}$ with non-zero elements from $i=[-68,68]$ and centered at $x=0$. The odd $W_{652}$ wavelet transform is applied to $\vec{G}^{(10)}$ to obtain the $G_{(m=1/2,k=8),i}$ coefficients, the even $W_{652}$ wavelet transform is applied to obtain the $G_{(m=1/2,k=7),i}$ coefficients. We obtain another gausslet denoted $\vec{G}^{(10,6)}$ also centered at $x=0$ with lower resolution by applying the $W_{652}$ scaling transform to the original $\mathcal{G}_{10}$. This process is repeated, by applying the odd, even, and scaling transforms $\vec{G}^{(10,6)}$, to to obtain $G_{(m=1/2,k=6),i}$, $G_{(m=1/2,k=5),i}$ and $\vec{G}^{(10,6,6)}$ respectively, which are now spaced $27\delta x$ apart. The resolution is further decreased using the same three transforms (odd, even and scaling) to obtain $G_{(m=1/2,k=4),i}$, $G_{(m=1/2,k=3),i}$ and $\vec{G}^{(10,6,6,6)}$ with spacing $81\delta x$. After the final set of transforms, we obtain $G_{(m=1/2,k=2),i}$,$G_{(m=1/2,k=1),i}$ and $\vec{G}^{(10,6,6,6,6)}$ with spacing $243{\scriptstyle\Delta}x$. Then $G_{m=0,k=0,i}$ is defined as $\vec{G}^{(10,6,6,6,6)}$. To obtain other $m$ values, one simply has to shift the coefficients by the appropriate number of positions for each resolution level. Namely $G_{(m,k),i}=G_{(m=1/2,k),i+(m-1/2)\times s}$ where $s=9$ for $k=7,8$, $s=27$ for $k=5,6$, $s=81$ for $k=3,4$ and $s=243$ for $k=1,2$. For $k=0$ $G_{(m,k),i}=G_{(0,0),i+m\times 243}$. The number of non-zero elements in each row of ${{{\bf G}}}$ is $252,706,2094,4528$ for each level respectively.

The calculation of the matrix elements is then simply a contraction of the full operator matrices as

[TABLE]

where ${{{\bf O}}}$ is one of ${{{\bf T^{(g)}}}}$ or ${{{\bf C^{(g)}}}}$ described by Eq. (70) and Eq. (71) respectively, and ${{{\bf G}}}$ has elements $G_{ni}$ as generated from wavelet transforms and translations. Clearly, very large matrices are needed in one dimension to calculate matrix elements in this fashion. However, the matrices are banded due to all matrix elements depending on the banded $S_{i^{\prime}i}^{(g)}$ such that the calculation and storage requirements only grow linearly with increasing number of wavelets.

The other issue is that a larger number of quadrature points are required to calculate accurate energy levels. With the Wilson basis, as few as seven Gauss-Jacobi quadrature points are needed for the quadrature error to be less than 0.001$a_{0}$. For the gausslet basis, $51$ Gauss-Legendre quadrature points are required to achieve the same precision. If a smaller $\delta x$ spacing is used than even more quadrature points are required. Also, Gauss-Jacobi quadrature does not work for gausslets, most likely due to the discontinuity in $V_{i^{\prime}i}^{(g)}$ causing problems here but are not important in the Wilson basis. Another issue is that the basis is not variational with respect to quadrature points. If fewer points than necessary are used, one can obtain energies that are lower in energy than the exact value.

The advantage of the gausslet/wavelet basis comes from trying to describe two different eigenfunctions simultaneously. As can be seen from Fig. 11, the accuracy of the 1S state is only marginally worse (for a given basis size) if both the 1S and 2S state are optimized simultaneously, as opposed to only optimizing the 1S state. To obtain chemical accuracy for only the 1S state requires 1366 basis functions, while obtaining chemical accuracy for both the 1S and 2S states requires 1717 basis functions.

This result is very different from what is found with the Wilson basis. To obtain chemical accuracy only optimizing the 1S state requires only 823 functions while obtaining chemical accuracy of both 1S and 2S states requires 2974 functions. Therefore, using gausslets should allow fewer necessary basis functions for describing multiple states at once. Symmetry adapted wavelets were used in the calculation along with Hermitian matrix-vector products.

IX.2 Augmenting the Wilson basis with Gaussian Slater type orbitals

Another possible improvement is the ability to combine the Wilson basis with well developed quantum chemistry basis sets. As an example, we will find the ground state of one-electron in a potential of a Li nuclei at the origin and two H nuclei at $\pm 3a_{0}$. The first thing to note is that the STO-$n$G is a basis in sum-of-products form such that

[TABLE]

Next, we can expand each of the $\exp\left[-(x-R_{i})^{2}\right]$ Gaussians in terms of the Wilson basis. We wish to describe the portion of phase-space that has contributions from all three nuclei with the Wilson basis set. Therefore, we partition phase-space such that the high momentum region ($k$ greater than cutoff $k_{c}$) will be described by the projection of the STO-$n$G basis into the $k>k_{c}$ Wilson basis functions. Thus, for $k>k_{c}$ the STO-$n$G basis is expanded as

[TABLE]

where $s^{(i)}_{m_{x}^{\prime},k_{x}^{\prime}}=\sum_{n_{x}=(m_{x},k_{x})}S^{-1/2}_{n^{\prime}_{x}n_{x}}\int_{-\infty}^{\infty}\exp\{-\alpha_{i}\left(x-R_{i}\right)^{2}\}d_{m_{x},k_{x}}\left(x\right)$ and $p_{m,k}(x)=\sum_{n_{x}}S^{-1/2}_{n_{x}n^{\prime}_{x}}d_{m_{x}^{\prime},k_{x}^{\prime}}\left(x\right)$. These states are then combined with the Wilson basis set for $k<k_{c}$.

These functions are orthogonal to the $k\leq k_{c}$ Wilson basis but are not orthogonal to each other, therefore orthogonalization needs to be performed in order to use this new basis iteratively. This is done using Gram-Schmidt starting with the smallest $\alpha_{i}$ and working upwards. The small values of $\alpha$ correspond to diffuse, delocalized orbitals while large values of $\alpha$ correspond to tight, localized orbitals. Thus, the larger values of $\alpha$ will contain higher momentum components. The Gram-Schmidt orthogonalization procedure preserves the ordering of the momentum when the input states are sorted by $\alpha_{i}$.

We can separate the sum of products and use individual components, $\sum_{m_{z},k_{z}>k_{c}}s_{m_{z},k_{z}}^{(i)}p_{m_{z},k_{z}}\left(z\right)$ representations as 1D basis functions. For this calculations, $n=6$ and we use the projected STO-$n$G basis functions only for the high momentum part of the Lithium nuclei. The two smallest values of $\alpha_{i}$ are diffuse and the part of phase-space these functions describe is described by the Wilson basis. This means that the basis used for the calculation is now

[TABLE]

where $i^{\prime}$ labels the Gram-Schmidt orthogonalized states. Applying this basis to the H-Li-H one-electron system with an even symmetrized basis increases the rate of convergence substantially. Using the basis of Eq. (77) required a basis size of 1493 to obtain the value of -5.1642 $E_{h}$ while using only the Wilson basis required 4250 basis functions. This comparison was made using ${\scriptstyle\Delta}x=0.5\sqrt{\pi}$ with 13 Gauss-Jacobi quadarture points and Hermitian matrix vector products.

The use of the STO-$n$G basis’ sum-of-products form presented here would most likely be even more useful if one represented the wavefunction as a sum-of-products basis. The reason is that the region of phase-space represented by the STO-$n$G basis only requires the addition of $n$ basis functions. Here, the $n$ basis functions had to be added in each dimension and then coupled into the rest of the basis so the savings was much smaller. That being said, the improvement in basis size is still substantial. The main downside is that calculating the necessary two-body Coulomb matrix elements for a fixed grid is no longer independent of the nuclear configuration. This is due to centering of the STO-$n$G basis functions on the nuclei.

IX.3 Other extensions

There are a few directions one could pursue to increase the efficiency of the Wilson basis: further exploits of the locality, the symmetry, and extending the extrapolation methods to estimate complete basis set limits.

First, it would be beneficial to take advantage of the fact that the individual matrices in the sum-of-products expansion of the Coulomb operator are banded in the sense that the off-diagonal elements decay exponentially. This is somewhat non-trivial to leverage as the bandedness depends on the quadrature point the matrix is evaluated at. In both Eq. (42) and Eq. (44), $p_{\alpha}$ closer to minus one is more diagonal in position space while $v$ closer to one is more diagonal in momentum space. Taking advantage of the bandedness would speed up the calculation and also assist with the troublesome memory usage of $L^{4}$.

Second, there room to further exploit the symmetries of the problem. As mentioned in Section VII, the translational symmetry in position space for the two-electron matrix elements has not been exploited here . It may also be useful to use a linear combination of Wilson basis functions that have arguments of the pairs of coordinates $(x_{1},x_{2}),(y_{1},y_{2})$ and $(z_{1},z_{2})$ instead of only $x_{i},y_{i}$ and $z_{i}$. First, symmetry adaptation of the basis would be more effective since it respects the symmetry of the fully correlated wavefunction. Second, as there is no obvious way to separate the $x_{1},x_{2}$ matrix of Eq. (44), the use of these basis functions would not decrease the speed of the calculation.

Third, expanding the basis using the technique of Section V.1 produced convergence of the calculated energies that can be fitted to a function with form $aM^{b}$. This should make it possible to extrapolate to the complete basis set limit although this was not done in the present study.

The use of phase-space localized (PSL) basis functions has garnered some attention to assist in the calculation of vibrationalHalverson and Poirier (2012); Shimshovitz and Tannor (2012b); Brown and Carrington (2016) energy levels. The motivating idea is that only a finite number of PSL basis functions would be needed in order to cover the localized region of phase-space in which the wave function occupies. For vibrational calculations, a comparison with PSL functions and the commonly used Gauss-Hermite functionsBrown and Carrington (2016) or sinc basisLarsson, Hartke, and Tannor (2016) suggests that this motivating idea has complications. It is important to describe the tail portion of the wavefunction that tunnels outside the classically allowed region of phase-space, but this tail region requires a much larger PSL basis. For vibrational problems, describing the tail region properly is difficult a priori as the potential becomes more complicated away from the minimum where a multi-dimensional Taylor expansion is accurate. This is not the case for electronic calculations where the potential can be expanded succinctly away from the nuclei using multipole expansions. Thus, in the future, it may prove beneficial to replace the tail regions $\left|m\right|>m_{c}$ with projections of the Wilson basis into standard basis sets.

X Conclusion

Is the Wilson basis competitive with commonly used basis sets at this time? The answer to this question is no. However, there are three advantages worth highlighting. First, this basis performs better in terms of computation and convergence than the previously investigated sinc basis functions. Second, for isolating a single state, the Wilson basis uses fewer functions than gausslets in the examples tested. Third, and most promising, we have shown how to combine the Wilson basis functions with the commonly used Slater type Gaussian basis functions.

Future areas of investigation for the Wilson basis include electronic systems with large applied magnetic fields, translationally invariant electronic systems, and use of the Wilson basis in quantum simulation algorithms on quantum computers. In order to model systems with an applied magnetic field, the extension to complex basis functions must be done and this is relatively straight-forward. Our preliminary computations suggests that the inclusion of a magnetic field does not greatly increase the number of basis functions needed for convergence although ${\scriptstyle\Delta}x$ must be adjusted.

In this paper, we have only looked at small molecular systems but there is an opportunity to study the electronic properties of solids with the Wilson basis functions. Since the Wilson basis functions are phase-space localized, the techniques used in this paper can equally be applied to calculations in momentum space.

A final application area of the Wilson basis function is quantum simulation algorithms for quantum computing. To utilize the power of quantum computation for electronic structure theory, it is crucial to map fermions to qubits optimally. Efficient mapping allows the number of overall quantum gates applied to be reduced, which is especially important given the current limitations of quantum computing hardware. While there are multiple fermion-to-qubits mappings Bravyi and Yu Kitaev (2000); Seeley, Richard, and Love (2012); Verstraete and Cirac (2005); Ball (2005); Bravyi and Yu Kitaev (2000); Havlíček, Troyer, and Whitfield (2017); Setia and Whitfield (2018), our most recent work has highlighted the potential of the Bravyi-Kitaev Super-Fast mapping Setia and Whitfield (2018) as well as its robustness against certain quantum noise processes Setia et al. (2018). This mapping, which has connections to lattice gauge theory Zohar and Cirac (2018), simulates the gauge fields rather than standard mappings that simulate the fermionic fields themselves. Because of the close connect to lattice gauge theories, our group has noticed strong dependence on the choice of basis set Chien et al. (2019). Future work will be in applying the Wilson functions described here to fermion encodings that require localized wave functions.

Acknowledgements

We gratefully thank the NSF for support under award number 1820747.