Adaptive algorithm for electronic structure calculations using reduction of Gaussian mixtures
Gregory Beylkin, Lucas Monzon, Xinshuo Yang
TL;DR
This paper introduces an adaptive algorithm for electronic structure calculations that reduces Gaussian mixtures efficiently, avoiding basis errors and requiring fewer terms than traditional multiresolution methods.
Contribution
The novel adaptive method selects the optimal Gaussian mixture terms during iteration, eliminating the need for multiresolution bases and improving efficiency.
Findings
Successfully applied to Hartree-Fock equations for diatomic molecules
Achieves comparable accuracy with fewer Gaussian terms
Demonstrates effective adaptivity without multiresolution bases
Abstract
We present a new adaptive method for electronic structure calculations based on novel fast algorithms for reduction of multivariate mixtures. In our calculations, spatial orbitals are maintained as Gaussian mixtures whose terms are selected in the process of solving equations. Using a fixed basis leads to the so-called "basis error" since orbitals may not lie entirely within the linear span of the basis. To avoid such an error, multiresolution bases are used in adaptive algorithms so that basis functions are selected from a fixed collection of functions, large enough as to approximate solutions within any user-selected accuracy. Our new method achieves adaptivity without using a multiresolution basis. Instead, as a part of an iteration to solve nonlinear equations, our algorithm selects the "best" subset of linearly independent terms of a Gaussian mixture from a collection that is…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3| 1 | ||
| 2 | ||
| 3 | ||
| 4 | ||
| 5 | ||
| 6 | ||
| 7 | ||
| 8 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Adaptive algorithm for electronic structure calculations using reduction
of Gaussian mixtures
Gregory Beylkin, Lucas Monzón and Xinshuo Yang
Department of Applied Mathematics
University of Colorado at Boulder
UCB 526
Boulder, CO 80309-0526
Abstract.
We present a new adaptive method for electronic structure calculations based on novel fast algorithms for reduction of multivariate mixtures. In our calculations, spatial orbitals are maintained as Gaussian mixtures whose terms are selected in the process of solving equations.
Using a fixed basis leads to the so-called “basis error” since orbitals may not lie entirely within the linear span of the basis. To avoid such an error, multiresolution bases are used in adaptive algorithms so that basis functions are selected from a fixed collection of functions, large enough to approximate solutions within any user-selected accuracy.
Our new method achieves adaptivity without using a multiresolution basis. Instead, as part of an iteration to solve nonlinear equations, our algorithm selects the “best” subset of linearly independent terms of a Gaussian mixture from a collection that is much larger than any possible basis since the locations and shapes of the Gaussian terms are not fixed in advance. Approximating an orbital within a given accuracy, our algorithm yields significantly fewer terms than methods using multiresolution bases.
We demonstrate our approach by solving the Hartree-Fock equations for two diatomic molecules, and LiH, matching the accuracy previously obtained using multiwavelet bases.
.
1. Introduction
We present a new adaptive method for electronic structure calculations based on algorithms for reduction of multiresolution multivariate mixtures [12]. While we represent solutions using a linear combination of Gaussians (which has a long history in quantum chemistry), these Gaussians are not selected in advance to be used as an approximate basis but are generated in the process of solving equations. Thus, we avoid the so-called “basis error” usually associated with methods using Gaussians. Our approach can also be characterized as a “gridless” or a “meshfree” method.
Using Gaussians to find solutions of quantum chemistry problems has its origins in seminal papers [13, 24, 29] and is motivated by the fact that integrals involving these functions can be evaluated efficiently. In these early quantum chemistry papers the authors proposed to use linear combinations of Gaussians whose exponents and coefficients were found (or optimized) via Newton’s method in order to capture the correct behavior near the nuclear cusps and the correct rate of decay. However, this approach proved unsustainable as problems became larger. Instead, the construction of a basis for spatial orbitals has been performed off-line and the resulting sets of functions were 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 a basis.
The use of these Gaussian bases have revolutionized computational quantum chemistry in spite of the absence of a systematic way for controlling accuracy or providing guaranteed error bounds. In fact selecting an appropriate basis set became an art form requiring insight into the underlying solution. However, once a basis set is selected, the accuracy of the solution obtained using such a basis is ultimately limited. Since the equations being solved, e.g. the Hartree-Fock equations or the Kohn–Sham equations of density functional theory (DFT) (see e.g. [26]), only provide approximate solutions, the limitation in accuracy makes it difficult to separate the impact of using approximate equations from the approximate methods for solving them. For this reason, the use of adaptive methods to avoid the loss of accuracy caused by basis sets is highly desirable in this field.
The advent of multiresolution analysis (see e.g. [14, 3]) laid a conceptual foundation for adaptive methods but it took some time before practical adaptive algorithms were developed using multiwavelets [1]; see [19, 31, 32]. A central element for the success of the multiresolution algorithms can be traced to the fact that physically significant (integral) operators arising in problems of quantum chemistry are naturally represented by radial kernels, which, in turn, can be accurately approximated by a linear combination of Gaussians. The key advantage of using such approximations is that they yield a separated representation in which operators are efficiently applied one direction at a time. Without a separated representation, multiresolution operators would be too expensive to be practical in dimensions three and higher.
Multiresolution methods systematically refine basis functions (i.e., numerical grids) in the vicinity of the cusp-type singularities while using relatively few basis functions elsewhere. These methods have proven successful in efficiently computing highly accurate solutions, achieving guaranteed error bounds and eliminating the basis error. This approach has been implemented at ORNL in the software package MADNESS ( Multiresolution ADaptive Numerical Environment for Scientific Simulation, see [18]), which is now considered the most accurate approach in this field [20]. Wavelets are also used in BigDFT [16] (cf. a mixed-basis method using plane waves and atom-centered radial polynomials [25] and interlocking multi-center grids [27], as examples).
However, since wavelets (or multiwavelets) do not resemble the spatial orbitals, a large number of basis functions is needed to represent solutions, e.g. an individual orbital may require basis functions in three dimensions. Moreover, such local refinement schemes do not take advantage of the essential simplicity of the spatial orbitals far from the nuclei and require boundary conditions to limit the computational domain. While adaptive multiresolution methods are sufficiently fast to be used within one-particle theories of quantum chemistry, an advancement towards two particle theories or solving Schrödinger’s equation had a limited success due to computational costs (see e.g. [7, 8]).
Hence it is of interest to develop new adaptive schemes and compare them with MADNESS. As it was demonstrated in [4], it is possible to iteratively solve equations of quantum chemistry using new algorithms based on nonlinear approximation of functions to reduce the number of terms in intermediate representations (without resorting to Newton’s method). The results of this approach are presented in [4] and are based on using Slater-type orbitals.
In this paper we present a new adaptive method that uses linear combinations of Gaussians to represent the solutions. We demonstrate the new approach by solving the Hartree-Fock equations for and LiH so that we can compare the resulting representations with those in [19, 4]. As in [19, 4], we formulate the problem using integral equations and use a convergent iteration to solve them. In contrast to the above approaches, we use linear combinations of Gaussians to accurately approximate not only operators and potentials but also the functions on which they operate. As a result, all integrals can be computed explicitly and exactly by simply updating the parameters of the Gaussians involved. The computational effort thus moves from that of approximating and computing integrals to that of maintaining a reasonable number of terms in intermediate representations of solutions within the iterative scheme. Using our approach, computing a single integral (e.g. involving a Green’s function, potential and a wave function) can easily yield a Gaussian mixture with terms, most of which are linearly dependent (within a user-selected accuracy). We use the reduction algorithm from [12, Algorithm 1] to reduce the number of terms after each operation by finding the “best” linearly independent terms, thus maintaining a reasonably small number of them. This algorithm has computational complexity , where is the original number of terms and is the number of so-called skeleton terms selected by the algorithm and is the cost of computing the inner product between Gaussians as a function of dimension (in this paper and is a small constant). There is no underlying grid to maintain (thus, “gridless” or “meshfree” method) and there is no need to impose boundary conditions to reduce the size of the computational domain as in [19].
Since we only present results using examples of two small molecules, it is natural to ask if our approach generalizes to large molecules. The obvious obstacle is the quadratic dependence of the reduction algorithm on the number of skeleton terms (clearly, the final number of terms to represent solutions for large molecules is expected to be large). The key to addressing this problem is to avoid the unnecessary application of the reduction algorithm to sets of functions that are obviously linearly independent. Algorithmically, it means splitting terms of a Gaussian mixture into groups that are likely to have linear dependence (and never perform a global reduction). This approach allows us to maintain a reasonable number of skeleton terms within each group. Moreover, the resulting reduction algorithm is trivially parallel since the reduction within each group is independent of other groups. We demonstrate subdivision into groups using both location and scale in our examples and show that the accuracy of the solution improves gradually due to a convergent iteration.
We also observe that the Hartree–Fock setup does not take into account the electron-electron cusps whereas it is well known that incorporating them significantly improves the accuracy of energy calculations. In our approach, it would be natural to introduce a linear combination of Gaussian geminals (originally proposed in the seminal papers [13, 29]). Using Gaussian geminals to improve accuracy has a long history in quantum chemistry (see [23, 2] on this topic and references therein). We note that incorporating geminal functional forms into the adaptive approach of this paper appears possible (but by no means trivial) since the reduction algorithm for multivariate mixtures discussed in [12] (and on which the approach of this paper is based) is applicable to much more general functions than those used in our two examples. We also note that a different approach to incorporate cusps and to address basis set error was recently proposed in [17] by using a range-separated DFT formulation. We plan to address adaptivity involving electron-electron interactions and to demonstrate our approach on large molecules elsewhere.
We start by briefly stating the Hartree–Fock equations in Section 2 and demonstrating that the functional form of our approximate solution, a Gaussian mixture, is maintained when solving the integral form of these equations via iteration (we note that the same reasoning applies to the Kohn-Sham equations). Then we briefly review approximations of operators and potentials that we use in Section 4. We then turn to the reduction algorithm in Section 3 and describe a subdivision scheme which is critical to the practical use of our approach. We present examples and comparisons in Section 4 and discuss the results in Section 5.
2. Solving the Hartree-Fock equations
We use an integral version of the Hartree–Fock equations (see e.g. [30]) and present these equations for the reader’s convenience in order to emphasize a combination of properties that make our approach possible. We note that our approach is equally applicable to the Kohn-Sham equations considered in [19]. Our method is based on computing orbitals as Gaussian mixtures and using the reduction algorithm in [12, Algorithm 1] to keep the number of terms under control.
Briefly, the occupied orbitals are the lowest eigenfunctions of the Hartree–Fock operator
[TABLE]
where , . The eigenvalues are referred to as orbital energies and are negative for the occupied orbitals. The external potential accounts for attraction of electrons to the nuclei at locations , ,
[TABLE]
where is the number of nuclei, is the charge of the th nucleus and is the standard Euclidean vector norm. The Coulomb operator describes the potential created by orbitals ,
[TABLE]
and the exchange operator is defined as
[TABLE]
where indicates the complex conjugate of the function . The orbitals are obtained solving the coupled equations (2.1) and the total energy, , is calculated by adding the nucleus–nucleus repulsion energies to the total electron energy
[TABLE]
In order to obtain integral equations for (2.1), we follow [21, 22, 19, 4] and use that the operator has Green’s function
[TABLE]
That is,
[TABLE]
Note that the kernel of in the equations above is
[TABLE]
and that (2.1) is equivalent to
[TABLE]
where . We solve (2.3) via an iteration in which are changing and approach . There are four observations (properties) that make our approach possible as follows.
- (1)
The iteration to solve (2.3) is convergent (see additional comments below). 2. (2)
All potentials and Green’s functions can be accurately and efficiently approximated by linear combinations of Gaussians. This property of potentials and Green’s functions is the foundation of the approach in MADNESS, see e.g. [19, 18]. 3. (3)
Since in our approach we represent orbitals via linear combinations of Gaussians and potentials and Green’s functions are approximated by linear combinations of Gaussians, all integrals are evaluated explicitly. 4. (4)
In order to keep the number of terms in the representations of the orbitals under control, the key tool in our approach is a reduction algorithm that selects a subset of the “best” linearly independent terms from a large number of terms which result from applying operators.
The integral form of the Kohn-Sham equations has the same properties (cf. [19], where properties 1 and 2 are essential) so that our approach can be used for these equations as well.
2.1. An iteration to solve integral
equations
We solve (2.3) via the following iteration. We initialize using a collection of Gaussians centered inside the convex polyhedron defined by nuclear centers (see Remark 4).
Then, at step , we first apply the operator to the current approximation of the orbitals and compute
[TABLE]
where
[TABLE]
[TABLE]
and
[TABLE]
(complex conjugation is removed noting that in our calculations orbitals are real). Next, we compute entries of the matrix
[TABLE]
and evaluate its eigenvalues, which we set to be approximations of the orbital energies , . We then compute
[TABLE]
and a new set of functions
[TABLE]
Finally, we orthogonalize (and normalize) the resulting functions , , yielding the updated orbitals , .
We observe that as long as the Green’s functions, the Coulomb potentials and the orbitals are represented as Gaussian mixtures, then each operation described above maintains this form, i.e. the result is again a Gaussian mixture. However the number of terms grows rapidly and we show how to control the number of terms in Section 3.
Remark 1*.*
This iteration is known to be convergent for computing bound states (see discussion in [19]) although we are not aware of a rigorous mathematical proof of this fact. An argument can be made that if the potential is in the so-called Rollnik class (see [28]) then, for any fixed , , the norms of relevant operators are less than . This makes this iteration a fixed point iteration. The nuclear Coulomb potentials in just miss being in the Rollnik class because of their slow decay away from the nuclei. However, any truncation of these potentials that is sufficient for computing the bound states (at an arbitrarily large distance from nuclei) puts them back into the Rollnik class. In this paper we simply accept the fact that the iteration is convergent. We do truncate the nuclear Coulomb potentials at infinity in our approximations (see below) as is done (explicitly or implicitly) in all numerical methods for computing bound states.
2.2. Accurate approximation
of Green’s functions and potentials
We briefly recall the key approximation results in [9, 10] for the power functions , , and the Green’s function for the non-oscillatory Helmholtz equations,
[TABLE]
where is the modified Bessel function of the second kind.
Theorem 2**.**
[10, Theorem 5]** For any , and , there exist a step size and integers and such that
[TABLE]
where
[TABLE]
and is any number .**
The error estimates are based on discretizing the integral
[TABLE]
using (an infinite) trapezoidal rule
[TABLE]
where the step size satisfies
[TABLE]
and is any user-selected accuracy. The proof in [10, Theorem 5] allows one to shift the grid by any , .
For a given accuracy , power and a range of values , the infinite sum (2.9) is then truncated due to the exponential or super-exponential decay of the integrand at , to yield a finite sum approximation in that range. The number of terms is estimated as
[TABLE]
This approximation provides an analytic construction of a multiresolution, separated representation for the Poisson kernel and the Coulomb potential in any dimensions. In practice many terms with small exponents in (2.8) can be combined using [9, Section 6] to reduce the total number of terms further. In our examples we use the following approximation of the Poisson kernel and the Coulomb potential via a linear combination of Gaussians:
[TABLE]
where , , and exponents and coefficients are given in Table 1.
A similar result holds for the Green’s function for the non-oscillatory Helmholtz equation,
[TABLE]
where the integrand has a super-exponential decay at . An accurate and efficient approximation is obtained by discretizing this integral representation. We use
[TABLE]
where the step size is selected to achieve the desired accuracy ,
[TABLE]
Note that this integral is discretized for a range of the parameter so that only the weights (coefficients) change as we change . In our examples we approximate the Green’s function (2.11) with
[TABLE]
resulting in for in the range .
3. Reduction algorithm with subdivision
of terms into groups
We approximate orbitals using Gaussian mixtures (cf. the multiwavelets in [19, 31, 32] and the Slater-type orbitals in [4]). While using a fixed basis set of Gaussians selected in advance has a long history in quantum chemistry, in contrast with previous methods we do not choose a set of Gaussians in advance (nor do we use a particular multiresolution basis of Gaussians although such an approach is possible using results in [11]). We allow the iterative algorithm for solving the equations in Section 2.1 to “select” the necessary basis functions. After each operation that yields a significant increase in the number of Gaussians, we prune the resulting Gaussian mixture by using the reduction algorithm in [12, Algorithm 1] that chooses the “best” linearly independent subset of terms, the so-called skeleton terms. By estimating the error of the solution, we terminate the iteration when the desired accuracy is achieved.
Conceptually our method is simpler and, also, less technical than either [19, 31, 32] or [4] and, as far as we know, is novel. Clearly, the key to enabling our approach is the reduction algorithm described and analyzed in [12]. This reduction algorithm has complexity , where is the initial number of terms, is the number of skeleton terms and is the cost of computing the inner product between the terms of a Gaussian mixture ( is a small constant in dimension ). If the number of skeleton terms, , is small, the computational cost of the reduction algorithm is completely satisfactory. However, we cannot assume that the number of skeleton terms (which are used to represent orbitals) is small for a large molecule.
In this paper we introduce an important modification of the reduction algorithm to avoid an excessive computational cost due to the quadratic dependence on the number of skeleton terms. The modification stems from a simple observation that, if two sets of functions are known to be linearly independent in advance, an attempt to find skeleton terms is wasteful. In many applied problems a multivariate mixture with a large number of terms can be split into subgroups of terms so that each subgroup reflects the local behavior of the function. For example, it is well understood that orbitals have cusps at nuclear centers but otherwise are smooth functions (in fact it is this structure of solutions that underpins the success of multiresolution methods such as [19, 31, 32]). Well-localized terms associated with cusps decay rapidly away from a nuclear center so that these terms will be linearly independent (or even nearly orthogonal) from similar terms at other nuclear centers. This suggests an acceleration technique for the reduction algorithm based on a hierarchical subdivision of the terms of the mixture by scale and location. Given such a subdivision, reduction is performed only within a local group. In such an approach the reduction of terms within each group is independent from other groups and, thus, the computation is trivially parallel.
In [12] we consider Gaussian terms in the form of -normalized multivariate normal distribution ,
[TABLE]
where the vector defines its location and the symmetric positive definite matrix controls it shape. In the setup of this paper, all matrices are diagonal and, therefore, functions represented as Gaussian mixtures are in separated form (see [5, 6]). In fact, in our examples it is sufficient to use rotationally symmetric Gaussians so that their shape is controlled by a single scalar parameter.
We subdivide terms by grouping them by their scale (shape) and location. We consider a single (global) group of flat Gaussians. i.e. Gaussians that have a significant support around more than one nuclear center and control behavior of orbitals far away from nuclear centers. The rest of the terms we split into groups according to their proximity to nuclear centers using nuclear centers as “seeds” in a Voronoi-type decomposition: all terms are split into groups by their proximity to the nuclei, i.e. a term located at belongs to a group associated with a given nucleus at location if it is the closest to it among all other nuclei, i.e. . Such subdivision ensures that the cost of reduction is proportional to the number of nuclei since the number of skeleton terms associated with each group is expected to be roughly the same given the structure of functions they represent. In our experiments we observed that the number of skeleton terms in these groups is still fairly large (e.g. several thousands terms) so further subdivision by scale and location is appropriate. Towards this end, we associate with each nucleus center and its group the distance and subdivide each group further by considering two parameters: (i) the radial intervals , that restrict the locations of the terms; (ii) the ranges of exponents , where is a parameter (see discussion below), identifying the shape of the terms. Within each radial interval , terms may have significantly different shapes determined by the shape matrices and, thus, cannot be linearly dependent. Effectively, it is a multiresolution subdivision strategy.
In the examples in Section 4, we initialize the spatial orbitals so that every term in the mixture has a diagonal shape matrix , where is the identity matrix. This forces all Gaussians involved in further computations to be rotationally symmetric, i.e. to have of the same structure. As already stated, we first identify a global group by considering all terms with , where is a parameter. Terms assigned in this group are significant in the vicinity of more than one nuclear center and capture the far-field behavior. The rest of the terms we divide into groups by their proximity to the nuclei and scale (shape). Within each group associated with nuclei , we identify a Gaussian term as a member of a subgroup, denoted by if
[TABLE]
where . The index controls the location and the index controls the scale (shape) of a term. The role of these parameters is similar to those labeling basis functions in Multiresolution Analysis (MRA).
Remark 3*.*
In general, subdivision into groups may differ in detail from the one we use for our examples in this paper. This depends on the type of functions we want to represent. In particular, the number of indices controlling location of terms (and, therefore, subgroups) can be larger than what we use in our examples (e.g. there could be three such indices).
We note that, for large and , the corresponding group consists of terms with a small essential support and located (relatively) far away from the nucleus center . Such terms (and, therefore, groups) are not needed to represent orbitals since orbitals have cusps only at the nuclei. To reduce the number of unnecessary groups, we do not maintain such groups.
In order to capture the cusp behavior, we introduce in the vicinity of nuclei an additional set of groups , so that if
[TABLE]
In our computation we ignore any terms that are not included in the groups or for all nucleus locations . Selecting and , the total number of groups (including the global group) is . It is worth noting that, when reducing a Gaussian mixture according to the group subdivision, some of the groups might be empty.
To estimate how subdivision into groups speeds up computations, let us consider a multivariate mixture with terms where the number of linearly independent terms is . Let us assume that we can subdivide these terms into groups with an equal number of terms (for simplicity of the estimate). We assume that the number of linearly independent terms in each group is then at most , where . The cost of performing the reduction algorithm for one group is then and, thus, the total cost to reduce all groups is . Depending on , we get a significant speed-up factor and note that the reduction of independent groups of terms is trivially parallel (a property that we did not use it in our computations). With the subdivision into groups described above, the number of groups is naturally proportional to the number of nuclear centers. Since the solution has a cusp at the nuclear center, we can consider to be a constant (i.e., each cusp locally requires roughly the same number of terms). Hence, the reduction algorithm with subdivision into groups is linear in the number of nuclear centers.
Remark 4*.*
By examining operators in (2.4)-(2.7), we observe that if we start an iteration with a Gaussian mixture with diagonal shape matrices , then all resulting Gaussian mixtures involved in the computation will have diagonal shape matrices as well. Moreover, location of the terms is restricted to a box defined by positions of nuclei. We observe that, in (2.4)-(2.7), convolutions and multiplications produce new Gaussians (with different shifts and/or shape matrices). While convolution does not change the center of the resulting Gaussian terms, multiplication of two Gaussians generates a Gaussian centered at a different location. Specifically, we have for the product of two normal distributions
[TABLE]
where
[TABLE]
If , , and , then we have
[TABLE]
Hence, all three coordinates of the new center are in between the corresponding coordinates of and . This implies that all Gaussian terms resulting from iteration (2.4)-(2.7) will have their centers inside a box defined by the minimum and the maximum of individual coordinates of nuclear centers in (2.2) (in fact, this result holds in any dimension). Moreover, if the Gaussians are rotationally invariant, then all terms will have their centers inside a convex polyhedron defined by the nuclear centers . We always initialize the iteration using Gaussians centered inside the convex polyhedron defined by the nuclear centers.
We note that the restricted location of terms of a Gaussian mixture in our approach cannot be achieved when using wavelet or multiwavelet bases, which force a computational box that is much larger than the box defined by nuclear centers.
4. Examples of computations for two small molecules
4.1. Helium Hydride,
First we solve the Hartree–Fock equation for
[TABLE]
with the potential
[TABLE]
where , , and .
As described in Section 2.1, our basic approach involves recasting (4.1) as an integral equation which we solve iteratively. The iteration (see below) is convergent and, while the rate of convergence is slower than quadratic, it is sufficiently fast so that only a few dozen iterations are required (see discussion in [19] and references therein). We represent the spatial orbital as a Gaussian mixture
[TABLE]
that is gradually constructed via iterations. In each step of the iterative solution, the number of terms in representing grows and we use [12, Algorithm 1] to control their number.
We approximate potentials and and the Poisson kernel using (2.10) and the Green’s function in (2.11) using (2.12).
For the reader’s convenience, we explicitly describe the iteration for solving (4.1) ,
[TABLE]
where denotes convolution.
To start the iteration, we use the single Gaussian
[TABLE]
We choose as the error tolerance for the reduction of terms in each group and terminate the iteration if . After iterations we obtain the orbital energy and the total energy . Using the MADNESS software [15] yields orbital energy and total energy so that, using these values as reference, our computations have absolute error for the orbital energy and for the total energy. As a result of the iteration, the total number of terms representing the orbital is . This number of terms is significantly smaller than the terms required when using multiwavelets [19, 31, 32] and a factor larger than that in [4], where the number of Slater-type orbitals is . However, our approach is less technical than that in [4] and, for this reason, it is easier to extend it to deal with large molecules. In Table 2, we show the number of terms in the global group , the number of non-empty groups and the maximum , minimum and average number of terms in these non-empty groups in and at the final step.
In Figure 4.1, we display the spatial orbital on the line connecting the two nuclei locations and .
4.2. Lithium Hydride, LiH
For our second example, we consider the Hartree-Fock equations for lithium hydride, LiH. The Hartree–Fock equations in this case are
[TABLE]
where
[TABLE]
[TABLE]
and
[TABLE]
For lithium hydride we have , , and . We approximate potentials and and the Poisson kernel using (2.10) and the Green’s function in (2.11) using (2.12).
To solve (4.3), we initialize orbitals as
[TABLE]
There are many possible initializations, e.g., we can use an approximation to and via Gaussian mixtures. Importantly, the initial approximation should be chosen so that the initial orbital energies are negative. After orthonormalizing the functions , we compute the initial energies as the eigenvalues of the matrix defined in (2.5). We then set using (2.6) and update the orbitals via (2.7).
Finally, when the desired accuracy is achieved, we compute the total energy as
[TABLE]
In this example, we choose an error tolerance of for the reduction of terms in each group. We terminate the computation if for . After iterations, we arrived at the orbital energies and computed with absolute errors of and , respectively, when compared with the reference energies and obtained using the MADNESS software. The computed total energy has an absolute error of compared with the value evaluated in MADNESS. The total number of terms to represent orbitals and is and . In Table 3, we show the number of terms in the global group , the number of non-empty groups and the maximum , minimum and average number of terms in the non-empty groups in the Gaussian mixture representations of , , and at the final step.
The resulting number of terms is again much smaller (about 100 times smaller) than that needed using MADNESS software and by a factor larger than that in [4], where the number of terms to represent the Slater-type orbitals is and .
In Figure 4.2 we display the spatial orbitals and on the line connecting the two nuclei locations and .
Remark 5*.*
We implemented our code in Fortran90 and compiled it with the Intel Fortran Compiler version 19.0.1.144. All computations were performed on a single core (without any parallelization) of an Intel i7-6700 CPU at 3.4 GHz on a 64-bit Linux workstation with 64 GB of RAM. Currently it takes about seconds to solve the Hartree–Fock equations for and seconds to solve the equations for LiH. We made no attempt to optimize our implementation. While the reduction algorithm that splits Gaussian mixtures into groups (see Section 3) is trivially parallel, we did not take advantage of this property in the current version of our code. Since, for small molecules, the main cost is in reduction, we expect a speedup factor of about , given the number of groups in Table 3. Such an acceleration factor should bring the timing close to that of MADNESS. We plan to do a careful speed comparison with existing methods separately.
Remark 6*.*
The reduction algorithm we used [12, Algorithm 1] is capable of achieving near full double precision accuracy if certain computations are carried out in quadruple precision. In addition to solving (4.1) and (4.3) using the reduction algorithm implemented in double precision and obtaining orbital energies , we continued iterations and switched to the reduction algorithm implemented in quadruple precision to improve the accuracy. We observed that, within several additional iterations, the orbital energies improved by gaining additional accurate digits. However, as pointed out in [12], the quadruple precision implementation is more than times slower than the double precision version (if the same accuracy is required). We do not include these results since, by modifying the reduction algorithm, it should be possible to avoid using quadruple precision. However, we note that, in most practical computations, an accuracy of is typically sufficient.
5. Conclusions and further
work
We presented a new adaptive algorithm for computing the electronic structure of molecules. Our approach is novel since we do not select a basis in advance and let the iteration for solving a system of integral equations identify the necessary functions to achieve a desired accuracy. The reduction algorithm developed in [12] selects the “best” subset of linearly independent terms after each operation that generates a large number of terms. These terms are Gaussians (identified by their parameters) and, since potentials and Green’s functions are also represented via Gaussian mixtures, all integrals are evaluated explicitly. Hence, in the algorithm, we only update the parameters of Gaussian mixtures.
We introduce in this paper a hierarchical subdivision of terms of Gaussian mixtures so that the computational cost of reduction is proportional to the number of nuclear centers. Although we did not implement a parallel version of our approach, it is clear that such an implementation will have a significant impact on the speed since reduction within each group of terms is independent of the other groups and computations are also mostly independent between orbitals.
We note that our method is less technical than that used in MADNESS (see [18]) and in [4] since our only tool is the reduction algorithm. We plan further work on our approach to verify its performance on large molecules.
Acknowledgement
We thank the anonymous reviewers for their comments and encouragement.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] B. Alpert. A class of bases in L 2 superscript 𝐿 2 {L}^{2} for the sparse representation of integral operators. SIAM J. Math. Anal , 24(1):246–262, 1993.
- 2[2] G. M. J. Barca and P-F. Loos. Three- and four-electron integrals involving gaussian geminals: fundamental integrals, upper bounds, and recurrence relations. The Journal of chemical physics , 147(2):024103, 2017.
- 3[3] G. Beylkin, R. Coifman, and V. Rokhlin. Fast wavelet transforms and numerical algorithms, I. Comm. Pure Appl. Math. , 44(2):141–183, 1991. Yale Univ. Technical Report YALEU/DCS/RR-696, August 1989.
- 4[4] G. Beylkin and T. S. Haut. Nonlinear approximations for electronic structure calculations. Proc. R. Soc. A , 469(20130231), 2013.
- 5[5] G. Beylkin and M. J. Mohlenkamp. Numerical operator calculus in higher dimensions. Proc. Natl. Acad. Sci. USA , 99(16):10246–10251, August 2002.
- 6[6] G. Beylkin and M. J. Mohlenkamp. Algorithms for numerical analysis in high dimensions. SIAM J. Sci. Comput. , 26(6):2133–2159, July 2005.
- 7[7] G. Beylkin, M. J. Mohlenkamp, and F. Pérez. Preliminary results on approximating a wavefunction as an unconstrained sum of Slater determinants. Proceedings in Applied Mathematics and Mechanics , 7(1):1010301–1010302, 2007.
- 8[8] G. Beylkin, M. J. Mohlenkamp, and F. Pérez. Approximating a wavefunction as an unconstrained sum of Slater determinants. Journal of Mathematical Physics , 49(3):032107, 2008.