## Abstract

We present a new method for electronic structure calculations based on novel algorithms for nonlinear approximations. We maintain a functional form for the spatial orbitals as a linear combination of products of decaying exponentials and spherical harmonics centred at the nuclear cusps. Although such representations bare some resemblance to the classical Slater-type orbitals, the complex-valued exponents in the representations are dynamically optimized via recently developed algorithms, yielding highly accurate solutions with guaranteed error bounds. These new algorithms make dynamic optimization an effective way to combine the efficiency of Slater-type orbitals with the adaptivity of modern multi-resolution methods. We develop numerical calculus suitable for electronic structure calculations. For any spatial orbital in this functional form, we represent its product with the Coulomb potential, its convolution with the Poisson kernel, etc., in the same functional form with optimized parameters. Algorithms for this purpose scale linearly in the number of nuclei. We compute electronic structure by casting the relevant equations in an integral form and solving for the spatial orbitals via iteration. As an example, for several diatomic molecules we solve the Hartree–Fock equations with speeds competitive to those of multi-resolution methods and achieve high accuracy using a small number of parameters.

## 1. Introduction

We present a new approach to electronic structure calculations based on recently developed algorithms for computing near-optimal exponential approximations of functions [1–5]. We maintain a functional form for the spatial orbitals consisting of linear combinations of products of decaying exponentials and spherical harmonics centred at the nuclear cusps. While such representations are similar to the classical Slater-type orbitals, in the course of computation, we optimize both the exponents and the coefficients in order to achieve an efficient representation of solutions and to obtain guaranteed error bounds. In this way, we combine the efficiency of traditional Slater-type representations with the adaptivity of current multi-resolution methods.

An approach of using nonlinear algorithms to find solutions of quantum chemistry problems has its origins in seminal papers [6–8]. In these papers, the authors used sums of Gaussians whose exponents and coefficients were optimized in a nonlinear manner in order to capture the correct behaviour near the nuclear cusps and the correct rate of decay. Similar approaches have been used (cf. [9,10]) to optimize the exponents in Slater-type approximations for diatomic molecules. Although such dynamic nonlinear optimization via traditional methods (e.g. Newton’s method) may produce a very efficient representation, it is practical only for very small molecules. Consequently, construction of spatial orbitals is traditionally performed offline, and the resulting sets of functions are then used as a fixed basis, leading to the so-called basis error if the actual solution is not well approximated within the linear span of such fixed basis.

While methods using bases of optimized sets of Gaussians have revolutionized computational quantum chemistry, they generally lack the ability to control the approximation accuracy in a systematic way and to guarantee the error bounds. Indeed, selecting a basis set is often an art form that requires insights into the underlying solution and, once it is selected, the accuracy of the solution obtained using such basis is ultimately limited. The limitation on accuracy of the traditional approach has spurred the development of multi-resolution methods (cf. [11–15]). Multi-resolution methods systematically refine numerical grids (or basis functions) in the vicinity of the cusp-type singularities while using a relatively few grid points (or basis functions) elsewhere. Multi-resolution methods have proved successful in efficiently computing highly accurate solutions and achieving guaranteed error bounds. However, these methods rely on approximating spatial orbitals using basis functions (e.g. piecewise polynomials in references [11–13]) which do not resemble the spatial orbitals that typically arise in quantum chemistry calculations. As a result, multi-resolution methods require many parameters to faithfully represent the cusps. Moreover, such local refinement schemes do not take advantage of the essential simplicity of the spatial orbitals far from the nuclei. While adaptive multi-resolution methods are sufficiently fast to be used within one-particle theories of quantum chemistry [11–15], their use in solving the multi-particle Schrödinger equation [16,17] is computationally costly.

In addition to multi-resolution methods based on multi-wavelets, let us also mention a mixed-basis method using plane waves and atom-centred radial polynomials (cf. [18]), as well as interlocking multi-centre grids (cf. [19]).

In the new approach, we use a functional form for the spatial orbitals that involves linear combinations of products of exponentials and spherical harmonics,
1.1that capture the nuclear cusps at **R**_{j} and the far-field behaviour. In contrast to standard approaches, the number and the values of the (potentially complex-valued) exponents and *α*_{lm,n} are not fixed in advance; instead, they are optimized throughout the course of the computation in order to achieve an efficient representation with a specified accuracy. For efficiency reasons, we also use an alternative representation for the spatial orbitals involving linear combinations of products of exponentials and spherical interpolating functions,
1.2where *K*_{n} are interpolating functions on the unit sphere associated with near-optimal quadratures (see references [20–22]). Specifically, the interpolating condition *K*_{n}(*ω*_{n′})=*δ*_{n,n′} holds at the spherical quadrature nodes which are invariant under the icosahedral group [20]. It turns out that the interpolating functions *K*_{n} can be expressed in terms of spherical harmonics (see §3*b*) and, therefore, these two representations are functionally equivalent. In §3*b*, we describe algorithms for converting between these two functional forms. We use one form or another depending on efficiency considerations as is typical in pseudo-spectral methods where the interpolating representation is more convenient for multiplying functions, whereas the spectral representation is more convenient for computing convolutions.

We develop a numerical calculus based on these functional representations. For example, for any two functions *ϕ*_{1}(**r**) and *ϕ*_{2}(**r**) that are already in this functional form, we develop algorithms to represent the product *ϕ*_{1}(**r**)*ϕ*_{2}(**r**) and the convolution *Δ*^{−1}*ϕ*_{1}(**r**) in the same functional form and with a small number of parameters. By casting the relevant electronic structure equations (e.g. the Hartree–Fock equations) in an integral form, we demonstrate that the functional forms in (1.1) and (1.2) can be used to solve for the ground states via iteration (using the framework developed in Harrison *et al*. [11,23–25]). The numerical calculus used within the iteration framework allows us to efficiently build up highly efficient representations for the solutions with guaranteed error bounds. It is also noteworthy that we compute the convolution (*μ*^{2}+*Δ*)^{−1}*ϕ*_{1}(**r**) in the spectral domain, which obviates the need to solve a large sparse matrix equation.

We consider the algorithms presented here to be of a preliminary nature, because we plan to explore a number of improvements and variants of this approach. However, our numerical experiments on diatomic molecules indicate that they are as efficient (e.g. for multiplying functions) or more efficient (e.g. for applying convolution kernels) than the corresponding operations in a multi-wavelet framework. The algorithms developed here also allow us to achieve high accuracy using a small number of parameters. In particular, fewer than 3000 complex-valued parameters are required for representing the spatial orbitals for molecules of helium hydride (HeH+) and of lithium hydride (LiH) with ≈5×10^{−7} absolute errors in the orbital energies; by contrast, standard multi-resolution methods require several orders of magnitude more parameters for a comparable accuracy.

In §2, we briefly review a method for computing near-optimal exponential representations [1–5] and some recently developed algorithms for interpolation and integration on the unit sphere using near-optimal quadratures (see references [20–22]). In §3, we introduce basic forms for representing solutions of quantum chemistry problems. We then develop, in §4, a numerical calculus based on the nonlinear representations (1.1) and (1.2) described in §3. We then use this numerical calculus in §5 to solve the Hartree–Fock equations for several diatomic molecules. Finally, we discuss some possible alternative formulations and directions for future research in §6.

## 2. Preliminaries

In this section, we collect some results required in the paper.

### (a) Near-optimal exponential representations

Given a decaying function *f*(*r*), *r*≥0, we use algorithms developed in references [1–5] to construct an exponential representation,
2.1
where *ϵ* is the desired approximation error and the exponents *α*_{m}, , and coefficients *a*_{m} are, in general, complex-valued. This representation uses a small number *M* of terms for the given approximation error *ϵ*, and results in highly efficient approximations. The approach to constructing such approximations has its origins in the so-called Adamjan–Arov–Kreĭn (AAK) theory for optimal rational approximations [26–28]. The version of the algorithm in this paper computes the exponents *α*_{m} and coefficients *a*_{m} from the 2*N*+1 equispaced samples *f*_{n}=*f*(*Rn*/(2*N*)), *n*=0,…,2*N*, where *R* and *N* are chosen so that *f*(*r*) is sufficiently sampled and |*f*(*r*)|<*ϵ*, *r*≥*R*.⇓

The basic steps of this construction are given below. In algorithm 1, *X*^{†} denotes the pseudo-inverse of the matrix *X* and *X*(*m*:*n*,:) denotes the submatrix consisting of rows *m* through *n* and *X*=(**x**_{1}⋯**x**_{M}) denotes the matrix consisting of the column vectors **x**_{j}.

We have implemented a fast singular value decomposition solver for step (1) using the randomized techniques in references [29–31] and the fact that Hankel matrices can be applied in operations via the fast Fourier transform. Our implementation reduces the overall cost of algorithm 1 to operations, where the (implicit) constant is small. Because typically *M*≪*N* (e.g. *M*=10), the cost of the algorithm is essentially linear in the number of samples *N*. In our experience, the choice of the singular value *σ*_{M}<*ϵ* in algorithm 1 always results in an error bound, as long as the function *f*(*r*) is smooth and decays rapidly (which is the case for the applications here). In fact, in the experiments discussed in §5, choosing a fixed value of *ϵ*=10^{−6} in algorithm 1 was sufficient to achieve ≈5×10^{−7} for the computed orbital energies in the Hartree–Fock examples. However, the approximation error should be checked *a posteriori* (using, for example, the already computed values *f*_{n}=*f*(*Rn*/(2*N*))), and a smaller singular value chosen if necessary. The connection between the accuracy *ϵ* and the *M*^{th} singular value *σ*_{M} in algorithm 1 is one of the key features of AAK theory [26–28].

We use algorithm 1 to compute exponential representations of radial functions *f*(∥**r**∥) that may or may not have a cusp at **r**=0. When the function *f*(∥**r**∥) does not have a cusp (e.g. *f*(∥**r**∥)=*e*^{−∥r∥2}), algorithm 1 produces the exponential representation to *f*(**r**) that is effectively smooth at **r**=0. Specifically, in the Taylor expansion for small ∥**r**∥,
2.2the odd moments (approximately) vanish to sufficiently high order. We note that this property naturally arises from application of algorithm 1 and is not imposed as an additional constraint.

### Remark 2.1

Although algorithm 1 suffices for computing exponential approximations when the approximation error *ϵ* is no smaller than ≈10^{−7} and *f*(*r*) is smooth and decays rapidly (which is sufficient for many electronic structure calculations), computing more accurate approximations using this algorithm may require quadruple precision.

For this reason, one of the possible alternative formulations described in §6 involves the so-called reduction algorithm [3]. Specifically, if the function *f*(*r*) is already a linear combination of *N* decaying exponentials, then the reduction algorithm constructs another representation of the same form, but with a smaller number *M*≪*N* of exponents. The basic idea behind this approach is that it is straightforward to construct a suboptimal representation of a given function as a sum of decaying exponentials (i.e. a representation that contains an excessive number of terms for a desired approximation error *ϵ*); the reduction algorithm may then be used to compute another exponential representation, but with a smaller number of terms. The reduction algorithm requires operations, and is therefore essentially linear in the number of suboptimal exponentials. Moreover, in contrast to algorithm 1, the reduction algorithm reliably yields (near) optimal representations with approximation error *ϵ* as small as 10^{−14}, and has high efficiency even when *f*(*r*) decays slowly or has 1/*r*-type singularities.

### (b) Integration and interpolation on the sphere

Let us briefly recall some results from Ahrens & Beylkin [20] on quadratures for efficient integration and interpolation on the sphere with nodes invariant under the icosahedral group.

Denoting the unit sphere in as , an orthonormal basis for is given by the spherical harmonics,
2.3where ** ω** is a unit vector, , the polar angle

*θ*∈[0,

*π*], the azimuthal angle

*ϕ*∈[0,2

*π*) and are the normalized associated Legendre functions, for

*m*≥0 and by for

*m*<0. The quadrature nodes constructed in Ahrens & Beylkin [20] are invariant under the icosahedral group and, thus, do not concentrate excessively near the poles (figure 1). These quadratures lead to a near-optimal integration of functions on the sphere, i.e. the number of nodes

*N*

_{L}required to integrate functions in the subspace of dimension (

*L*+1)

^{2}, where 2.4is the subspace of spherical harmonics of maximum degree

*L*, satisfy (

*L*+1)

^{2}/(3

*N*

_{L})≈1. Following Ahrens & Beylkin [20], let us consider the reproducing kernel for , 2.5with 2.6where

*P*

_{l}are the Legendre polynomials and is the Jacobi polynomial. Given a function and

*N*≥2

*L*of its samples at the quadrature nodes , we arrive at an interpolation formula for the sphere (an analogue of Lagrange interpolation), 2.7where Similar to the polynomial Lagrange interpolation, (2.7) is convenient for performing pointwise operations such as multiplication. By contrast, the representation using spherical harmonics is convenient for performing the convolutions with the (bound state) Helmholtz and Poisson kernels that arise in quantum chemistry.

We note that if a function of three variables in spherical coordinates for any fixed *r*≥0, then this function can be represented by *N* functions (*N*≥2*L*) of one variable, .

We also make use of the addition theorem [32], which states that for two unit vectors, , 2.8and, thus, the representation (2.7) can also be expressed in terms of spherical harmonics.

## 3. Basic forms for representing functions

### (a) Basic form of the representation

In our electronic structure calculations, we represent functions in the form of
3.1where *f*_{c}(**r**) captures the behaviour near the cusps and *f*_{s}(**r**) accounts for the remaining smooth part. We note that Pahl & Handy [18] also use a decomposition into smooth and cusp terms (although the functional forms, and associated algorithms, are different).

Specifically, we represent the cusp part *f*_{c}(**r**) in the form
3.2where **R**_{j}, *j*=1,…,*J* denote the positions of the nuclei. Here, the radial component satisfies
3.3and denotes the interpolating function on the sphere associated with *N*_{c} spherical quadrature nodes and interpolation order *L*_{c} (see §2*b* for details). For the sake of simplicity, we assume that the number of quadrature nodes, *N*_{c}, does not depend on the singularity location **R**_{j} (in the examples in §5, *N*_{c}=12). In contrast to representations based on Slater-type orbitals or Gaussian-type orbitals, the number and the values of (complex-valued) exponents in the radial representation (3.3) are *not* fixed in advance. Indeed, as described in §4, the exponents (and the coefficients) are determined dynamically throughout the course of the computation to achieve a desired level of accuracy while using a near optimally small number of terms. We also note that the *J* nuclei-centred terms that comprise the cusp part *f*_{c}(**r**) may overlap with one another (and, in fact, do overlap in the examples in §5).

Similarly, we represent the smooth part *f*_{s}(**r**) in the form
3.4where denotes the interpolating function on the sphere associated with *N*_{s} spherical quadrature nodes and interpolation order *L*_{s}. The exponents and coefficients are again determined dynamically using algorithm 1. We note that because the form of the function *f*_{s}(**r**) efficiently captures the far-field behaviour of solutions to the Hartree–Fock equation, we expect that *N*_{s} can also be taken reasonably small (in the examples in §5, *L*_{s}=11 and *N*_{s}=192 is large enough to achieve ≈5×10^{−7} errors in the Hartree–Fock orbital energies). In addition, by construction the exponential representations in (3.4) will have vanishing moments of high order, and *f*_{s}(**r**) will be effectively smooth at **r**=**0** (see the discussion leading to equation (2.2)). Finally, we note that, from the interpolation property ,
This simple observation leads to a fast evaluation of *f*_{s}(**r**) on spherical grids , and is crucial for the efficiency of the algorithms in §4.

As described in §3*b*, the cusp part *f*_{c}(**r**) may be equivalently expressed in terms of spherical harmonics,
3.5where the radial component is again represented in terms of exponentials,
In an analogous way, the smooth part *f*_{s}(**r**) can also be expressed in terms of spherical harmonics,
3.6The representations (3.2) and (3.4) based on spherical interpolating functions are convenient for algorithms such as multiplication (see §4*a,b*), whereas the forms (3.5) and (3.6) are convenient for applying convolution operators (see §4*c*).

### (b) Converting between interpolating and spherical representations

We make use of two basic forms for representing *f*_{c}(**r**) and *f*_{s}(**r**) and now show how to convert between these two representations. It suffices to consider
3.7and
3.8We make use of the fact that the interpolating function can be written in the form
3.9where and denote the *nth* node and weight associated with the *N*-point quadrature rule (see §2*b*).

In order to convert from (3.7) to (3.8), we have
where the last equality uses formula (3.9). Thus, to convert from a spherical harmonic representation to an interpolating representation, it suffices to use algorithm 1 and the above formula for evaluating *f*_{lm}(*r*).

Similarly, to convert from (3.8) to (3.7), we use
Hence, it suffices to use algorithm 1 and the above formula for evaluating *f*_{n}(*r*).

### Remark 3.1

Algorithm 1 for representing *f*_{lm}(*r*) as a sum of near-optimal exponentials requires computing a spherical harmonic transform, for each radial grid point *r*=*r*_{j}. This can, in principle, be accomplished using a fast algorithm. While we plan to address this issue elsewhere, in the examples in §5, the subspace dimension *L* and the number *N* of spherical quadrature nodes (necessary to achieve high accuracy) are small, e.g. *L*=11 and *N*=192, and the direct evaluation is sufficiently fast. We also note that the spherical transforms (one for each radial grid point) can be computed in a trivially parallel manner.

## 4. Algorithms for quantum chemistry applications

We now discuss algorithms for performing electronic structure calculations using the representations from §3*a*. For the applications that we consider, it is sufficient to develop algorithms for computing representations of *f*_{1}(**r**)*f*_{2}(**r**), *f*_{1}(**r**)+*f*_{2}(**r**), *V* (**r**)*f*_{1}(**r**), (−*Δ*+*μ*^{2})^{−1}*f*_{1}(**r**) (*μ*≥0), and *Δf*_{1}(**r**) in the form discussed in §3*a*. Here, *V* (**r**) denotes the Coulomb potential,
4.1We develop algorithms for the non-trivial operations, *f*_{1}(**r**)*f*_{2}(**r**), *V* (**r**)*f*_{1}(**r**) and (−*Δ*+*μ*^{2})^{−1}*f*_{1}(**r**), and *Δf*_{1}(**r**).

### (a) Multiplication of representations

Let us consider two functions *f*(**r**)=*f*_{c}(**r**)+*f*_{s}(**r**) and *g*(**r**)=*g*_{c}(**r**)+*g*_{s}(**r**) in the form (3.1) in §3*a*. We present an algorithm to construct the same type of representation for the product,
where *ϵ* is the desired approximation error.

The basic idea behind the algorithm is that, in a neighbourhood of the singularity at **R**_{j}, the product *f*(**r**)*g*(**r**) has the form
where is a polynomial and *h*_{Lc}(*r*** ω**) is continuous. Therefore, the difference
has

*L*

_{c}+1 derivatives at

**r**=

**R**

_{j}, and can be efficiently represented in the form (3.6); we also observe that the far-field behaviour of

*f*(

**r**)

*g*(

**r**) is efficiently captured by (3.6). Noting the equivalence of the spherical harmonic and interpolating representations (3.2) and (3.5), it is sufficient to construct functions

*h*

_{j}(

**r**),

*j*=1,…,

*J*, such that for

**r**close to

**R**

_{j}, It then holds that the function

*h*

_{s}(

**r**), is effectively smooth (neglecting terms of size ), and can be efficiently represented in the form 4.2In fact, for

**r**in a neighbourhood of a possible cusp location

**R**

_{j}, Because the functions

*h*

_{k}(

**r**),

*k*≠

*j*are smooth at

**R**

_{j}, so is

*h*

_{s}(

**r**) (neglecting terms of size ). We observe that choosing

*L*

_{c}=2 and

*N*

_{c}=12 yields an accuracy of ≈10

^{−7}for examples in §5, and requires about 200 (complex-valued) parameters for each cusp.

In order to construct functions *h*_{j}(**r**), we compute, via algorithm 1, exponential representations for each spherical quadrature node , *n*=1,2,…,*N*_{c},
4.3Here, we can choose the cut-off *R* to be small, for example (only local behaviour matters when subtracting off the singularity). Using interpolation on the sphere (see §2*b*), we then have that for *r* small enough and all ,
4.4

The computation of *h*_{s}(**r**) is similar. In particular, for each spherical quadrature node , *n*=1,2,…,*N*_{s}, we construct via algorithm 1 exponential representations
4.5It then holds that, for all *r*≥0 and ,
where we assume that the number *N*_{s} of nodes on the sphere is chosen large enough to sufficiently sample the smooth part. In using algorithm 1, we choose the cut-off *R*, so that the product is smaller than the desired approximation error *ϵ*.

Let us now comment on the efficiency of this approach. First, constructing the cusp part, *h*_{c}(**r**), of the product *f*(**r**)*g*(**r**) requires sampling the functions and on an equispaced grid for each spherical quadrature node . Because the number *N*_{c} of spherical quadrature nodes is small (e.g. *N*_{c}=12), constructing these samples is inexpensive. In addition, the exponential representations (4.3) need only be accurate for 0≤*r*≤*R*, where *R* is small (e.g. ) and, hence, the number of required radial samples is also small. Once these samples are computed, applying algorithm 1 is inexpensive, because its cost is essentially linear in the number of samples.

Similarly, constructing the smooth part, *h*_{s}(**r**), of *f*(**r**)*g*(**r**) requires sampling the functions and on an equispaced grid for each spherical quadrature node . Here, the number, *N*_{s}, of spherical quadrature nodes is typically a factor of 10 or so larger than *N*_{c} (e.g. *N*_{s}=192 in the experiments described in §5), and it may first appear that evaluating these samples would be more costly. However, because the smooth terms *f*_{s}(**r**), *g*_{s}(**r**) and *h*_{s}(**r**) of the multiplicands and the product share the same quadrature nodes , sampling these terms requires no interpolation. In particular, recall from §3*a* that *f*(**r**) is of the form *f*(**r**)=*f*_{c}(**r**)+*f*_{s}(**r**). Because *f*_{c}(**r**) involves only *N*_{c} terms, where *N*_{c} is small, evaluating is inexpensive. In addition, although *f*_{s}(**r**) nominally involves a sum of *N*_{s} terms, we can simply use the interpolation relationship,
so that *f*_{s}(**r**) is also inexpensive to sample.

In addition to the above considerations, we also note that the computation of the exponential representations in (4.4) and (4.5) can be performed in a trivially parallel manner for each spherical quadrature node and .

To illustrate the basic idea of algorithm, we consider the function
where **R**_{1}=(0,0,−.7) and **R**_{1}=(0,0,.7). We use a subspace dimension *L*_{c}=2 and *N*_{c}=12 spherical quadrature nodes for the cusp part, and an error tolerance of *ϵ*=10^{−6} in algorithm 1. For the smooth part, we use *L*_{s}=11 and *N*_{s}=192. In figure 2*a*, we display the function *f*(**r**)*g*(**r**) on the line **r**=(0,0,*z*) connecting the two cusp locations **R**_{1} and **R**_{1}. In figure 2*b*, we display the function *f*(**r**)*g*(**r**)−*h*_{1}(**r**)−*h*_{2}(**r**) after the cusp part is subtracted. Note that we allow the two cusp terms *h*_{1}(**r**) and *h*_{2}(**r**) to overlap in order to avoid (artificially) sharp gradients in the representation.

### (b) Multiplication with the potential

Given the function *f*(**r**)=*f*_{c}(**r**)+*f*_{s}(**r**), we now discuss how to construct a representation
where *V* (**r**) denotes the Coulomb potential (4.1) and the functions *v*_{j}(**r**) decay exponentially fast. The algorithm is nearly identical to the multiplication algorithm in §4*a*. The main difference is that the cusp part *g*_{c}(**r**) is represented in the form
where the radial component is given by
4.6The smooth part *g*_{s}(**r**) is again represented in the form
The advantage of maintaining this intermediate form is twofold. First, by incorporating the 1/*r*-type singularities explicitly in the representation of *g*_{c}(**r**), only a small number of parameters need to be maintained for an accurate approximation of *V* (**r**)*f*(**r**). Second, multiplication by the potential *V* (**r**) is always followed by an application of the Green function, (−*Δ*+*μ*^{2})^{−1}(*V* (**r**)*f*(**r**)), and the resulting (less singular) function can again be efficiently represented in the basic form discussed in §3*a* (see §4*c* for additional details). We note that, as it was observed in the algorithms [11–13], multiplication by the potential typically creates a large number of fine scales which are then immediately reduced by convolution, with some loss of the overall efficiency. By maintaining the above intermediate form, we avoid this issue.

Similar to §4*a*, the basic idea is that, in a neighbourhood of the singularity **R**_{j}, the product *V* (**r**)*f*(**r**) has the form
where is a polynomial and *h*_{L+1}(*r*** ω**) is continuous. Therefore, the difference
has

*L*

_{c}derivatives at

**r**=

**R**

_{j}, and can be efficiently represented in the form (3.6). Thus, to construct the representation (4.6), we proceed as in §4

*a*and first construct functions

*g*

_{j}(

**r**),

*j*=1,…,

*J*, such that 4.7Once

*g*

_{j}(

**r**)

*j*=1,…,

*J*are computed, we construct the smooth part as Putting all this together, we have

We now describe the computation of *g*_{j}(**r**) in more detail. For each spherical quadrature node , *n*=1,2,…,*N*_{c}, we construct, via algorithm 1, exponential representations,
where
Then, extending *g*_{j}(**r**) by interpolation,
we arrive at the representation (4.6). The computation of the smooth part *g*_{s}(**r**) proceeds in the same way as the computation of the smooth part of *h*_{s}(**r**) of the multiplicand *f*(**r**)*g*(**r**) in §4*a*.

### (c) Convolution with the bound-state Helmholtz and Poisson kernels

We first discuss how to evaluate the convolution operator,
4.8to obtain the representation in the form described in §3*a*. We assume that *f*(**r**) is already given in this form and, for now, *μ*>0. As is typical in pseudo-spectral methods, it is more efficient to apply this operator in the Fourier domain—in this case, the space of spherical harmonics. Thus, we first convert *f*(**r**) from an interpolating representation to a spherical harmonic representation, as described in §3*b*.

By linearity and translation invariance of (−*Δ*+*μ*^{2})^{−1}, it also suffices to consider *f*(**r**) in the form . Now, by expanding the kernel in (4.8) via spherical harmonics, it can be shown that
where *F*_{l}(*r*) is given by
4.9and *i*_{l}(*z*) and *k*_{l}(*z*) are defined in terms of the modified Bessel functions
Thus, computing a representation of (−*Δ*+*μ*^{2})^{−1}*f*(**r**) simply involves using algorithm 1 to construct an exponential representation of *F*_{l}(*r*). We also note that when *f*(**r**) is of the intermediate form (4.6), the only difference is that we also have integrals of the form
for which we also need to construct exponential representations.

Recall that algorithm 1 requires sampling *F*_{l}(*ρ*) on an equispaced grid *r*_{n}=*Rn*/(2*N*), *m*=0,…,2*N*−1; here *R* is chosen large enough that |*F*_{l}(*R*)|≤*ϵ*, where *ϵ* is the desired approximation error. To do so efficiently, we use the recursion
The integral between *r*_{n} and *r*_{n+1} can be accurately evaluated using a small number of quadrature points because *i*_{l}(*μρ*) does not oscillate for *μ*>0 (e.g. we use five quadrature nodes in our examples in §5). Once the samples *I*_{1}(*r*_{n}) are computed, samples of the first integral in (4.9) may be readily obtained. We note that the values of *i*_{l}(*μρ*) at the quadrature points may be computed using an interpolation table constructed off-line. In addition, by assumption the function *f*_{0}(*ρ*) is represented with a near optimally small number of exponents, and is thus also inexpensive to sample. The second integral defining *F*_{l}(*r*) may be efficiently sampled in a similar manner.

Computing a representation of *Δ*^{−1}*f* (i.e. in the case *μ*=0) proceeds in a similar manner. In particular, we have
where
4.10One technical difference is that *F*_{l}(*r*) decays like *A*_{l}*r*^{−l−1}, where *A*_{l} depends on the exponents and coefficients of *f*_{0}(*r*); thus, for small *l*, the direct application of algorithm 1 would require a prohibitively large number of samples for a representation on the entire half line *r*≥0. However, in the quantum chemistry applications considered here, we need to only represent products of the form *g*(*Δ*^{−1}*f*), where the function *g*(**r**) (and, thus, the overall product) decays exponentially fast. Therefore, it suffices to construct an exponential approximation to *F*_{l}(*r*) only within the numerical support of *g*(**r**), which is the approach taken in the examples in §5.

### Remark 4.1

We can alternatively use the methods in references [2,3] to construct an efficient exponential representation of *F*_{l}(*r*) on the entire half line *r*≥0. In fact, the slowly decaying asymptotic part *A*_{l}*r*^{−l−1} can be represented with a small number of decaying exponentials via the discretization of an appropriate quadrature formula [2]. Once this exponential representation for *A*_{l}*r*^{−l−1} is available, an exponential representation for the rapidly decaying function *F*_{l}(*r*)−*A*_{l}*r*^{−l−1} can be constructed using algorithm 1.

## 5. Example of solving Hartree–Fock equations for diatomic molecules

We now use the representations and algorithms in §§3*a* and 4 to solve the Hartree–Fock equations for several diatomic molecules.

### Example 5.1

As our first example of using the representations and algorithms discussed above, we solve the Hartree–Fock equation,
5.1with the potential
As in references [11–13], our basic approach involves recasting (5.1) as an integral equation that we solve via iteration. However, in contrast to references [11–13], we represent *ϕ*(**r**) as described in §3*a* and apply algorithms described in §4. Because the spatial orbital *ϕ*(**r**) has cusp-like singularities at the nuclei **r**=**R**_{1} and **r**=**R**_{2}, the Hartree–Fock equation provides a useful accuracy test for our approach.

Let us now describe the solution method in greater detail. We write (5.1) as
where *μ*^{2}=−2*E* and
We solve this via the following iteration:
The inner products,
are computed by representing *f*(**r**)*g*(**r**) in the basic form (3.5) and (3.6) using the algorithm in §4*a*, and computing the resulting integrals analytically (a similar comment applies to computing 〈*f*,*V* *g*〉). We let the above iteration run until the computed correction to the orbital energy *E* is less than the desired accuracy.

We solve equation (5.1) with *Z*_{1}=−1, *Z*_{2}=−2 and **R**_{1}=(0,0,−0.7) and **R**_{2}=(0,0,0.7), corresponding to the HeH+. We use *N*_{c}=12 spherical quadrature nodes for the cusp part, *N*_{s}=192 nodes for the smooth part and a tolerance of *ϵ*=10^{−7} in algorithm 1. We note that the number of radial samples required when using algorithm 1 is typically a factor of ≈10 more than the number of resulting exponential/coefficient pairs (and varies depending on the spherical quadrature direction).

We verify that, for this choice of parameters, the orbital energy *E*=−1.660544 is computed with six significant digits, and with 4.4×10^{−7} absolute error. The comparison is made using the MADNESS software [33] yielding the orbital energy *E*=−1.66054378. The total number of exponent/coefficient pairs for the representation of the orbital *ϕ*(**r**) (in the form (3.6)) is 637.

### Remark 5.2

We have implemented example 5.1 in Fortran 90 and made a preliminary speed comparison with the MADNESS code. MADNESS uses a similar iteration scheme, but is based on an adaptive basis representation via multi-wavelets [11]. In our comparison, we choose parameters in both codes to achieve relative error ≈10^{−8} for the orbital energy *E*. Comparing timings for this example, the approach in this paper is about 2.5 times faster than the MADNESS code. We made no attempt to optimize our implementation, in part, because the algorithms presented here are likely to be improved upon as we expect some modifications (outlined in §6) to speed up our initial implementation by an additional factor. We plan to do a careful speed comparison in a separate publication.

### Example 5.3

For our second example, we consider the Hartree–Fock equation for LiH. We have
5.2where ,
The iteration is the same as in the first example, with two modifications. First, the spatial orbitals *ϕ*_{1}(**r**) and *ϕ*_{2}(**r**) are orthogonalized after each iteration. Second, the orbital energies *E*_{1} and *E*_{2} are updated after each iteration by solving for the eigenvalues of the 2×2 matrix . In order to compute the inner products of the form , we use a straightforward modification of the algorithms in §4*a,b* to represent of *ψΔϕ* in the basic functional form of §3*a*. As in §4*b*, the 1/*r*^{2}-type singularities at **r**=**R**_{1} and **r**=**R**_{2} are carried explicitly. Once the representation for *ψΔϕ* is constructed, the integral is evaluated analytically.

In this example, we again use *N*_{c}=12 spherical quadrature nodes for the cusp part and *N*_{s}=192 nodes for the smooth part. We also take and , and again choose a tolerance of *ϵ*=10^{−7} in algorithm 1.

The computed orbital energies *E*_{1}=2.4517624 and *E*_{2}=0.2978234 agree to six significant digits with the values *E*_{1}=2.451763 and *E*_{2}=0.297823 computed by the MADNESS software [33], and have the absolute errors of 5.2×10^{−7} and 4.1×10^{−7}. The computed total energy *E*_{tot}=−7.986937 also agrees to six significant digits with the value *E*_{tot}=7.9869364 computed in MADNESS, and has an absolute error of 7.7×10^{−7}. The total number of exponent/coefficient pairs for the orbitals *ϕ*_{1}(**r**) and *ϕ*_{2}(**r**) (in the form (3.6)) is 1282 and 1327.

Figure 3 displays the spatial orbital *ψ*(**r**)≡*ϕ*_{2}(**r**), the cusp part *ψ*_{c}(**r**) and the smooth part *ψ*_{s}(**r**), on the line **r**=(*x*,0,0) connecting the two nuclei locations **R**_{1} and **R**_{2} (plots for the other spatial orbitals for LiH and for HeH+ are similar).

## 6. Alternative formulations

We are currently considering several possible modifications of our approach.

— It may be possible to modify our approach so that, instead of using algorithm 1, we apply the reduction algorithm [3]. This change will likely result in a faster algorithm, and is currently being explored.

— It may be desirable to use only the cusp part (i.e. (3.2) and (3.5)) in the functional representations; this can be achieved with minor modifications to the algorithms in §4 by using a partition of unity (see also reference [19]), where the function

*χ*_{j}(**r**) vanishes at the singularity locations**R**_{k}≠**R**_{j}. This approach may also result in more efficient representations and algorithms.— It may be possible to speed up the algorithms by performing the nonlinear approximations on a sparse subset of the spherical quadrature nodes (and treat the resulting collection of exponentials as a fixed basis for the remaining spherical quadrature nodes). Effectively, this implies a multi-resolution approach for computing on the sphere, a development that has many additional applications.

— Instead of using decaying exponentials in (1.1) and (1.2), we can use decaying Gaussians with complex-valued exponents. Although such a switch will increase the number of terms in the orbital representations, working with Gaussians should simplify the computations in a number of ways. We plan to examine this trade-off.

## Funding statement

This research was partially supported by NSF grant no. DMS-100995 and DOE/ORNL grant no. 4000038129.

## Acknowledgements

We thank Dr Robert Harrison for sharing with us his insights into computational quantum chemistry, for his many helpful suggestions and comments, and for his assistance in using the MADNESS software. We also thank Dr Martin Mohlenkamp for his careful reading of the manuscript, and for the help in comparing our method with that based on multiwavelets.

- Received April 12, 2013.
- Accepted July 25, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.