In recent years, interest in the random-phase approximation (RPA) has grown rapidly. At the same time, tensor hypercontraction has emerged as an intriguing method to reduce the computational cost of electronic structure algorithms. In this paper, we combine the particle-particle random phase approximation with tensor hypercontraction to produce the tensor-hypercontracted particle-particle RPA (THC-ppRPA) algorithm. Unlike previous implementations of ppRPA which scale as O(r6), the THC-ppRPA algorithm scales asymptotically as only O(r4), albeit with a much larger prefactor than the traditional algorithm. We apply THC-ppRPA to several model systems and show that it yields the same results as traditional ppRPA to within mH accuracy. Our method opens the door to the development of post-Kohn Sham functionals based on ppRPA without the excessive asymptotic cost of traditional ppRPA implementations.
I. INTRODUCTION
The random phase approximation (RPA) was originally derived in the context of a homogeneous electron gas over 60 years ago.1,2 Although there are multiple forms of the RPA, research has largely focused on the particle-hole RPA (which is often referred to in the literature simply as the “RPA” without qualification). Despite some initial success, its application to quantum chemistry was largely overshadowed by the successful development of coupled-cluster theory in the 1970s and 1980s. However, the last 15 years have seen a renewal of interest in the RPA.3–6 The RPA is attractive in part because it can capture long-range dispersive interactions,7,8 has no static correlation error,9–12 and does not experience divergence for materials with small band-gaps.5 One of the more interesting theoretical developments to emerge from the study of the RPA was the realization that the phRPA could be rewritten in terms of coupled-cluster theory.13 This correspondence was used by several authors13,14 to develop efficient algorithms for the solution of the phRPA equations.
More recently, our group has examined the particle-particle RPA (ppRPA).15–19 Initial results indicate that, like the phRPA, the ppRPA has several interesting properties as a density functional, and is more accurate than the phRPA for a variety of chemical problems.15,18 The ppRPA successfully captures the known fractional-spin and fractional charge constraints on the exact functional, as seen in its correct dissociation of
Like the phRPA, the ppRPA can be rewritten as a coupled-cluster-like equation, in this case one which only includes so-called “ladder” terms.16,20 One shortcoming of the ppRPA is that the particular form taken by its ladder-coupled-cluster-like equations make it resistant to the techniques used to achieve an O(r4) scaling for the phRPA, where r is the size of the basis set. Whether the ppRPA is written in its matrix form, in its coupled-cluster form, or in its linear, double-integral form,15,19 all methods used for its solution thus far require O(r6) operations, provided that the number of holes and particles both scale as O(r). In this paper, we will show that a recently developed form of tensor compression can reduce this scaling from O(r6) to O(r4) with minimal reduction of accuracy.
Tensor compression decomposes high-rank tensors, which are abundant in electronic structure theory, into lower rank tensors. One well-known example of tensor compression is the resolution-of-identity (RI) technique for the compression of the 2-electron repulsion integral (ERI) tensor
RI methods such as RI-SVS or RI-V express this rank-4 tensor as the product of rank-3 and rank-2 tensors21–25
If the number of auxiliary functions PR = r(r + 1)/2, then Eq. (2) is an equality; if PH < r(r + 1)/2, then the equality in Eq. (2) is only approximate. However, the locality of the atomic orbitals makes it possible to obtain a very good approximation with only a small number of auxiliary functions, usually PR ≈ 4r.
Tensor hypercontraction goes one step further than RI techniques in that it expresses high-rank tensors entirely in terms of rank-2 tensors.26–28 For example, the THC expression for the ERI tensor is
As in the case of RI techniques, the use of PH = r(r + 1)/2 auxiliary functions will guarantee an exact representation of the ERI tensor. In practice, it has been shown that only PH = O(r) auxiliary functions are needed to obtain a good approximation to the ERI tensor, due again to the locality of the atomic orbitals.26
The THC of the ERIs can be obtained in two steps. First, we select a desired atomic basis (e.g., cc-pVDZ) and apply the standard RI-V method to decompose the ERIs in terms of PR auxiliary functions in the corresponding auxiliary basis (e.g., cc-pVDZ-RI) from the EMSL basis set library,29 yielding a decomposition of the form in Eq. (2). This step can be performed in O(r4) operations. We then perform a nonlinear fitting procedure to approximate the rank-3 tensor
Tensor hypercontraction can be applied to other objects in electronic structure theory besides the ERI tensor, especially to excitation tensors
The remainder of our paper will be organized as follows: in Sec. II, we will describe the THC-ppRPA algorithm. In Sec. III, we will present the results of applying the THC-ppRPA algorithm to several sets of small molecules, hydrogen chains, and alkanes, comparing our results to the exact ppRPA results obtained through traditional methods. We will show that only a linear number of auxiliary functions are needed to obtain good accuracy and that, when a linear number of auxiliary functions is used, the algorithm does indeed scale as O(r4). In Sec. IV, we will present our conclusions and areas for future work.
II. ALGORITHM
A. Traditional ppRPA
There are several equivalent formulations of the ppRPA equations,15 but we will use the formulation in terms of a Ricatti equation, as described in Refs. 14,20. After some rearrangement, we find that the solution of the ppRPA equations requires us to solve
Here, εa and εi are the orbital energies of the unoccupied and occupied orbitals a and i in the reference state. Note that this reference state can be chosen to be either the Hartree-Fock state or a Kohn-Sham reference state using some particular DFT functional. The use of Kohn-Sham reference systems can be justified based on the formulation of time-dependent DFT on the pairing field.33 Although this approach is not common practice in coupled cluster theory, DFT references are often used in Green's function methods, including the RPA. When we solve Eq. (4) for the excitation tensor
When spin is taken into account, the RPA equations for closed-shell singlet systems separate into two different channels: a same-spin αα channel and an opposite-spin αβ channel. The ββ channel is identical to the αα channel since the α and β spins are equivalent for singlets. Furthermore, the solution ααT of the same-spin equations is simply the antisymmetrized version of the solution αβT of the opposite-spin equations such that the total correlation energy can be written as
where the second term captures the effects of exchange and could be dropped to implement an exchange-free DFT functional. Because of this simplification, we will first calculate the solution αβT of the opposite-spin channel and then use this solution to calculate the energy in both the opposite-spin and same-spin channels. Although there is no reason that our method cannot be applied to non-singlet states, we will assume for the remainder of the paper that we are dealing with closed-shell singlets to avoid the costs and complications of higher spin states.
In the traditional implementation of the algorithm, Eq. (4) can be solved iteratively in much the same way that coupled-cluster equations are solved. Assuming that
Set
$^{(0)}T^{ab}_{ij} = 0$, n = 1.- Calculate(7)\begin{eqnarray}{}^{(n)}{T}^{ab}_{ij} &=& -\frac{1}{\epsilon _a+\epsilon _b-\epsilon _i-\epsilon _j}\big(\epsilon ^{ab}_{ij}+{}^{(n-1)}{T}^{ab}_{cd}\epsilon ^{cd}_{ij}\nonumber\\&&+ \epsilon ^{ab}_{kl} {}^{(n-1)}{T}^{kl}_{ij}+{}^{(n-1)}{T}^{ab}_{kl}\epsilon ^{kl}_{cd} {}^{(n-1)}{T}^{cd}_{ij}\big).\end{eqnarray}
If ‖(n)T − (n − 1)T‖ is sufficiently small or some maximum number of iterations has been reached, then stop. Otherwise, set n → n + 1 and return to step 2.
Although this method will work, it will unavoidably scale as O(r6) because step 2 requires the matrix multiplication of ε and T, both of which are O(r2) × O(r2) matrices. To achieve a speed-up using THC methodology, we need to avoid explicit evaluation or multiplication of the matrices ε and T. Fortunately, it is possible to adapt the iterative methodology outlined above to create an efficient algorithm, which we describe in Sec. II B.
B. THC methodology
To achieve a speed-up for our algorithm, both ε and T must be written in THC form. There exist numerous methodologies for writing the ERI tensor
Next, we have to represent the excitation tensor
In Eq. (8), the matrix zPQ is symmetric such that the symmetry (a, i) ↔ (b, j) is preserved. We also recall that the antisymmetry of a ↔ b and i ↔ j need not be taken into account since we are solving only the opposite-spin channel; the correlation energy of the same-spin channel will be included at the end of the calculation when we evaluate Eq. (6).
As we said, in order to achieve an improvement in algorithmic efficiency, we must avoid any explicit evaluation of
where
But if we remember that both
To make our iterative procedure more stable, we add a damping term to Eq. (9) which tends to lead to smoother convergence. Our final cost function is
This damping approach is common in Hartree-Fock calculations and adds minimal computational cost to our algorithm. We find that a selection of p = 0.5 is effective at damping oscillations that can occur during iterative optimization.
The cost of our algorithm comes from evaluating the value and gradient of J during our minimization procedure. To see that this cost scales as O(r4), we can evaluate the various terms in J to see that they require the summation over at most four indices. For instance, one of the terms in our cost function will be
where we have defined
The sum in Eq. (12) can be evaluated in O(r4) time, provided that the number of auxiliary functions scales as PA = O(r). In contrast, evaluating Eq. (12) without using tensor hypercontraction requires the summation over the six indices a, b, c, d, i, j and O(r6) time. If we similarly expand the other terms in Eq. (11), we find that J and its gradients dJ/dx, dJ/dy, dJ/dz can all be evaluated in O(r4) time, provided that both PA and PH scale as O(r).
Based on our discussion above, the THC-ppRPA algorithm can be outlined as follows:
III. RESULTS
We applied our THC-ppRPA algorithm to several molecules in a variety of basis sets. First, we considered the calculated correlation energy of six small molecules in the cc-pVDZ basis as a function of PA, the size of the auxiliary basis set. Table I shows the results. As the number of auxiliary functions is increased, the correlation energy approaches the traditional ppRPA result. By the time the number of auxiliary functions PA is less than twice the size of the basis set size r, the correlation energy is within 1 mH of the traditional result, far less than the error of the traditional ppRPA algorithm itself.
The THC-ppRPA correlation energy in mH for six small molecules in the cc-pVDZ basis as a function of the number of auxiliary functions. The final column lists the result of traditional ppRPA. In all cases, only PA = 40 auxiliary functions were required for the THC-ppRPA result to converge to within a few mH of the traditional ppRPA correlation energy.
. | . | . | . | . | . | Exact . |
---|---|---|---|---|---|---|
Molecule . | PA = 20 . | PA = 30 . | PA = 40 . | PA = 50 . | PA = 60 . | ppRPA . |
CH2(r = 25)) | 75.85 | 79.64 | 79.90 | 79.92 | 79.98 | 79.83 |
CO (r = 30)) | 198.80 | 211.68 | 214.78 | 214.18 | 214.59 | 214.54 |
H2O (r = 25) | 148.79 | 155.33 | 154.83 | 154.97 | 154.97 | 154.76 |
HCN (r = 35) | 178.75 | 201.77 | 207.88 | 208.99 | 209.03 | 208.68 |
N2 (r = 30) | 209.47 | 224.03 | 227.84 | 228.22 | 227.89 | 227.69 |
NH3(r = 30) | 128.80 | 139.21 | 140.60 | 140.48 | 140.39 | 140.21 |
. | . | . | . | . | . | Exact . |
---|---|---|---|---|---|---|
Molecule . | PA = 20 . | PA = 30 . | PA = 40 . | PA = 50 . | PA = 60 . | ppRPA . |
CH2(r = 25)) | 75.85 | 79.64 | 79.90 | 79.92 | 79.98 | 79.83 |
CO (r = 30)) | 198.80 | 211.68 | 214.78 | 214.18 | 214.59 | 214.54 |
H2O (r = 25) | 148.79 | 155.33 | 154.83 | 154.97 | 154.97 | 154.76 |
HCN (r = 35) | 178.75 | 201.77 | 207.88 | 208.99 | 209.03 | 208.68 |
N2 (r = 30) | 209.47 | 224.03 | 227.84 | 228.22 | 227.89 | 227.69 |
NH3(r = 30) | 128.80 | 139.21 | 140.60 | 140.48 | 140.39 | 140.21 |
We can also examine convergence properties of the THC-ppRPA algorithm with respect to the number of iterations. Figure 1 shows the behavior of the THC-ppRPA algorithm with PA = 60 when applied to H2O in the cc-pVDZ basis. Because neither the traditional ppRPA algorithm nor the THC-ppRPA algorithm is variational, the THC-ppRPA result will not approach the ppRPA result monotonically. However, the correlation energy does eventually converge to the traditional result. Figure 1 also shows the benefit of including an “anti-damping” term in our cost function, which smoothes out oscillations in the convergence.
The convergence of the THC-ppRPA algorithm with and without a damping term for H2O in the cc-pVDZ basis.
The convergence of the THC-ppRPA algorithm with and without a damping term for H2O in the cc-pVDZ basis.
We observe the same results when we apply our algorithm to the same six molecules in the larger cc-pVTZ basis set. Table II shows the correlation energies in this larger basis. Once again, auxiliary sets of approximately PA = 2r functions are sufficient to recover a good approximation to the traditional correlation energy.
The THC-ppRPA correlation energy in mH for six small molecules in the larger cc-pVTZ basis as a function of the number of auxiliary functions.
. | . | . | . | . | . | Exact . |
---|---|---|---|---|---|---|
Molecule . | PA = 30 . | PA = 50 . | PA = 70 . | PA = 90 . | PA = 110 . | ppRPA . |
CH2(r = 65)) | 93.03 | 104.91 | 107.87 | 108.27 | 108.34 | 108.49 |
CO (r = 70)) | 210.79 | 256.53 | 278.91 | 282.59 | 283.68 | 283.03 |
H2O (r = 65) | 176.11 | 204.42 | 207.83 | 209.89 | 210.15 | 209.36 |
HCN (r = 85) | 176.35 | 234.78 | 258.33 | 270.04 | 272.56 | 272.79 |
N2 (r = 70) | 232.95 | 276.52 | 292.39 | 295.83 | 297.08 | 296.01 |
NH3(r = 80) | 154.20 | 176.52 | 172.35 | 185.14 | 186.59 | 186.09 |
. | . | . | . | . | . | Exact . |
---|---|---|---|---|---|---|
Molecule . | PA = 30 . | PA = 50 . | PA = 70 . | PA = 90 . | PA = 110 . | ppRPA . |
CH2(r = 65)) | 93.03 | 104.91 | 107.87 | 108.27 | 108.34 | 108.49 |
CO (r = 70)) | 210.79 | 256.53 | 278.91 | 282.59 | 283.68 | 283.03 |
H2O (r = 65) | 176.11 | 204.42 | 207.83 | 209.89 | 210.15 | 209.36 |
HCN (r = 85) | 176.35 | 234.78 | 258.33 | 270.04 | 272.56 | 272.79 |
N2 (r = 70) | 232.95 | 276.52 | 292.39 | 295.83 | 297.08 | 296.01 |
NH3(r = 80) | 154.20 | 176.52 | 172.35 | 185.14 | 186.59 | 186.09 |
These examples are interesting and suggest that the number of auxiliary functions needed to obtain good accuracy is far less than the theoretical limit of PA = r(r + 1)/2. However, they are not clear evidence that the THC-ppRPA algorithm provides an asymptotic speed-up over the traditional approach because the number of electrons in a particular molecule is constant as the basis set size is increased. To obtain evidence of genuine asymptotic scaling improvement, we first examine linear hydrogen chains in the cc-pVDZ basis. In these systems, we can tabulate the fraction of the traditional ppRPA correlation energy recovered by the THC-ppRPA as a function of the number of auxiliary functions. Table III shows our results. Here, we can clearly see that the number of auxiliary functions PA needed to obtain any particular fraction of the correlation energy (say, 90%) scales linearly with r rather than quadratically. The same result is seen in alkane chains, as shown in Table IV. Once again, the number of auxiliary functions needed scales as r. This result implies that we are achieving an asymptotically significant compression of the excitation tensor
The percentage of correlation energy recovered by THC-ppRPA algorithm relative to the traditional ppRPA algorithm for hydrogen chains of length N in the cc-pVDZ basis. This table shows that the number of auxiliary functions required to achieve some constant fraction ppRPA energy scales roughly linearly with the size of the system.
Molecule . | PA = N . | PA = 2N . | PA = 3N . | PA = 4N . | PA = 5N . |
---|---|---|---|---|---|
H2 | 57.8% | 90.5% | 94.6% | 99.8% | 100.1% |
H4 | 50.8% | 79.2% | 93.7% | 98.1% | 99.8% |
H6 | 41.2% | 75.0% | 92.7% | 98.2% | 99.8% |
H8 | 43.3% | 73.6% | 92.9% | 97.6% | 99.4% |
H10 | 43.4% | 74.0% | 92.3% | 97.0% | 99.2% |
H12 | 40.4% | 71.0% | 91.0% | 96.8% | 99.2% |
H14 | 46.8% | 74.8% | 91.2% | 95.9% | 98.9% |
H16 | 44.2% | 69.7% | 90.5% | 95.9% | 98.8% |
Molecule . | PA = N . | PA = 2N . | PA = 3N . | PA = 4N . | PA = 5N . |
---|---|---|---|---|---|
H2 | 57.8% | 90.5% | 94.6% | 99.8% | 100.1% |
H4 | 50.8% | 79.2% | 93.7% | 98.1% | 99.8% |
H6 | 41.2% | 75.0% | 92.7% | 98.2% | 99.8% |
H8 | 43.3% | 73.6% | 92.9% | 97.6% | 99.4% |
H10 | 43.4% | 74.0% | 92.3% | 97.0% | 99.2% |
H12 | 40.4% | 71.0% | 91.0% | 96.8% | 99.2% |
H14 | 46.8% | 74.8% | 91.2% | 95.9% | 98.9% |
H16 | 44.2% | 69.7% | 90.5% | 95.9% | 98.8% |
The percentage of correlation energy recovered by THC-ppRPA algorithm relative to the traditional ppRPA algorithm for alkanes with N carbon atoms in the cc-pVDZ basis. Again, the number of auxiliary functions required to achieve some constant fraction ppRPA energy scales roughly linearly with the size of the system.
Molecule . | PA = 20N . | PA = 40N . | PA = 60N . | PA = 80N . |
---|---|---|---|---|
CH4 | 87.2% | 99.7% | 100.1% | 100.2% |
C2H6 | 89.2% | 99.6% | 99.9% | 100.2% |
C3H8 | 90.2% | 99.7% | 99.8% | 100.0% |
C4H10 | 91.1% | 99.4% | 100.0% | 100.1% |
Molecule . | PA = 20N . | PA = 40N . | PA = 60N . | PA = 80N . |
---|---|---|---|---|
CH4 | 87.2% | 99.7% | 100.1% | 100.2% |
C2H6 | 89.2% | 99.6% | 99.9% | 100.2% |
C3H8 | 90.2% | 99.7% | 99.8% | 100.0% |
C4H10 | 91.1% | 99.4% | 100.0% | 100.1% |
This linear scaling of PA with r leads to an overall asymptotic reduction in the complexity of the algorithm, since the evaluation of the cost function and its derivatives require summation over at most 4 indices. We have also verified this scaling numerically. Figure 2 shows the time required for a single evaluation of the cost function and gradient as a function of the size of the hydrogen chain or alkane molecule along with the linear best-fit of the data. The O(r4) scaling of the algorithm is clearly visible in the case of the hydrogen chain. The slight deviation of the alkane running time scaling from the O(r4) ideal is likely due to the fact the number of electrons, holes, and basis functions do not grow at a rate exactly proportional to the length of the carbon chain. For instance, if the running time for alkanes were plotted against the number of electrons, a linear best fit would yield a scaling of r4.24.
The computational cost in ms of a single evaluation of the cost function and gradient as a function of the size of the basis set, for hydrogen chains and alkanes. Hydrogen chains show clear O(r4) scaling. Alkanes show slight deviation, which is primarily due to the fact that the number of electrons and holes are not quite directly proportional to the number of basis functions.
The computational cost in ms of a single evaluation of the cost function and gradient as a function of the size of the basis set, for hydrogen chains and alkanes. Hydrogen chains show clear O(r4) scaling. Alkanes show slight deviation, which is primarily due to the fact that the number of electrons and holes are not quite directly proportional to the number of basis functions.
IV. DISCUSSIONS AND CONCLUSIONS
In this article, we have introduced the THC-ppRPA algorithm for improving the asymptotic efficiency of the traditional ppRPA algorithm. Our method relies on the tensor decomposition of both the ERI tensor 2ε and the excitation tensor 2T. Results here and elsewhere confirm that the number of auxiliary functions needed to accurately approximate both of these objects scales as O(r) rather than the theoretical limit of O(r2), where r is the total number of basis functions in our system. This reduction of order in the number of required auxiliary functions leads to a reduced asymptotic scaling; whereas the traditional ppRPA algorithm scales as O(r6), our THC-ppRPA algorithm scales as O(r4).
Two important insights were necessary for the formulation of our algorithm. First, the ppRPA equation had to be rewritten as a coupled-cluster-like equation. We implemented several other variations of THC-ppRPA that relied on double-integral and matrix formulations of the ppRPA equation, but these were unsuccessful in providing a stable, efficient algorithm. Only the coupled-cluster formulation allowed us to solve the ppRPA equations efficiently within a THC expansion. Second, the coupled-cluster equation had to be recast as an iterative, scalar optimization problem. Solving the original Ricatti equation in Eq. (4) directly would not lead to a speed-up because it requires us to obtain the full, rank-4 tensor 2T from its THC components, an expansion that requires O(r6) time. However, we could have converted Eq. (4) into a non-iterative scalar optimization problem by minimizing the cost function
Although this approach would have allowed us to avoid explicitly expanding T, notice that it requires us to calculate the term
At this point, it is helpful to put our work in the context of the broader application of THC to coupled-cluster theory. It was shown in Ref. 28 that CCSD can theoretically be implemented in O(r4) operations using tensor hypercontraction. Because algorithms like CISD, QCISD,28 and CC232 can all be obtained through a truncation of terms in the CCSD equations, the result in Ref. 28 proves that the scaling of these other algorithms can also be reduced through the application of THC. What is more, the recent discovery that phRPA and ppRPA can both be written as coupled-cluster-like equations16,20 showed that THC-RPA can also be viewed—in principle—as a subset of THC-CCSD. However, in practice, many of these THC methods must be treated separately from THC-CCSD because they will display quite different scaling behavior than the THC-CCSD algorithm and will make use of different solution techniques.
A good example of the differences that can arise even between closely related THC algorithms is provided by the contrast between phRPA and ppRPA. Both phRPA and ppRPA are—in principle—special cases of the CCSD equations. Both phRPA and ppRPA can also be straightforwardly implemented using exact matrix diagonalization that requires O(r6) operations.5 However, it was shown that the scaling of phRPA could be reduced to O(r4log r) and later to O(r3) by applying RI36 and THC techniques,14 respectively. Despite the close connection of the phRPA and ppRPA algorithms, the ppRPA equations cannot be solved as efficiently as the pHRPA equations due to the different ways in which the indices are entangled in the two different methods. For example, the most costly term in the Ricatti equations for both pHRPA and ppRPA involves the product T · ε · T (see the last term in Eq. (4). In ppRPA, this term can be written as the THC expansion,
Notice that there is no way to evaluate the quantity inside the brackets in fewer than O(r5) operations. In contrast, the corresponding term in the phRPA equations can be written as
In this case, the quantity inside the brackets can be evaluated in only O(r3) operations through a clever ordering of summations. This marked difference in complexity is what enables THC-phRPA to be implemented with O(r3) scaling while the best known scaling for THC-ppRPA is the O(r4) scaling demonstrated in this paper. This example shows that, despite their similarities, the phRPA and ppRPA equations cannot be solved using the same techniques with the same reduction in cost. The THC-ppRPA algorithm should be developed and studied separately from other THC algorithms, as we have done in this article.
One important question that arises from our work is when the crossover between traditional ppRPA and THC-ppRPA occurs. For what system size r does THC-ppRPA become more efficient than traditional ppRPA? Here, we are hesitant to offer a definitive number for multiple reasons. First, there are many different equivalent approaches for solving the traditional ppRPA equations (matrix diagonalization, double integral formulation, Ricatti formulation) and each will require different computational resources. For example, our own experience has shown that the solution of the double integral formulation (see Ref. 15) can be tens to hundreds of times more costly than the solution of the matrix formulation. Because the application of the ppRPA method to molecules is relatively new, it is not yet clear which methodology should be viewed as the “standard” one or how it can be implemented most efficiently. The crossover point between traditional ppRPA and THC-ppRPA will therefore depend on which “traditional” ppRPA method is implemented. Second, the cost of the THC-ppRPA algorithm will depend on the number of auxiliary functions PA and PH required for an accurate approximation of the excitation tensor and the ERI tensor. Because this number will be strongly system-dependent, the crossover point will likewise be strongly system-dependent. A related concern is that the level of accuracy deemed acceptable will determine the number of auxiliary functions and the degree of convergence required, both of which will affect the THC-ppRPA algorithm's cost. For all these reasons, while we believe that the crossover point is likely to be large (r ≈ 500) due to the high prefactor of THC-ppRPA, declaring a single crossover point would depend on too many factors to be of much use. And certainly, future implementations of THC-ppRPA can work on reducing the crossover point without knowing the precise location of the crossover point for a particular system.
Two major areas for future work exist. The most important area for improvement is the large computational prefactor required by our algorithm. While we have shown that our asymptotic scaling is O(r4), the use of THC introduces a substantial prefactor that, for small values of r, can overwhelm any improvement due to the lower-order scaling. This result is not surprising, given that the same problem is faced by other THC algorithms. Nonetheless, it represents a very important challenge whose solution could significantly improve the practicality of THC-ppRPA. Currently, a single evaluation of the THC-ppRPA cost function and its gradient for a problem like H2O in the cc-pVDZ basis is not significantly more expensive than its solution by traditional ppRPA. However, THC-ppRPA requires many evaluations of the cost function and its gradient, first for the optimization of the cost function J during a single iteration and then for the repeated iterations required to solve the coupled-cluster equation. The overall performance of the algorithm could then be improved either by reducing the number of optimization steps or the number of iterations required for convergence.
A second area of investigation involves the use of THC-ppRPA to correct DFT functionals. While long-range dispersive forces are not captured by many popular DFT functionals, they are captured by RPA.5,15,18 The ppRPA also accurately describes the dissociation of
Speaking broadly, the applicability of THC to many different electronic structure algorithms is encouraging. It suggests that many if not all electronic structure algorithms can be improved by the application of tensor hypercontraction. Nonetheless, because each algorithm is distinct, it is worthwhile to explicitly develop THC versions of different algorithms to develop a larger toolkit for applying THC to conventional algorithms. For instance, the strategy for applying THC to CISD is different than for applying it to MP2. It is our hope that the extension of THC to the ppRPA algorithm will facilitate further application of THC to areas which have not yet been considered.
ACKNOWLEDGMENTS
N.S. would like to acknowledge support from the UNC EFRC: Center for Solar Fuels, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0001011. W.Y. acknowledges support from the National Science Foundation under award CHE-1362927.