What are the computational limits of explicitly correlated methods?

Question

Expanding the $N$-electron wavefunction in terms of Slater determinants (as in CI and CC theory) could lead to a very slow convergence to its basis set limit because such expansions can't give an accurate representation of the electron-electron cusps. As noted by Slater (and Hylleraas), going beyond determinants, including the interelectronic distance explicitly, can improve this.

This idea was exploited by Hylleraas and he wrote the wavefunction as $$ \Psi=e^{-\zeta(r_1+r_2)}\sum_{i=1}^{N}c_i(r_1+r_2)^{l_i}(r_1-r_2)^{2m_i}(r_{12})^{n_i}, $$ where $r_i$ is the distance of electron $i$ from the nucleus and $r_{12}$ is the distance between two electrons. For the helium, he found that a 308-term expansion gives an accuracy within $10^{-9}$ Hartree (!). However, in Jensen's Introduction to Computational Chemistry (2007 edition) it says that Hylleraas type wavefunctions are impractical for more than $3$ or $4$ electrons. However, there have been a lot of development in computational methods in the last $13$ years (improvement in parallel processing, GPUs, etc..), so, this can't be that bad.

Furthermore, in 1960, Boys and Singer found that functions that are product of Gaussian orbital and factors of the type $e^{-a r_{ij}^2}$ generate relatively simple integrals to calculate. But, in Piela's Ideas of Quantum Chemistry (2006) it says that the area of application of this method is relatively small due to computational reasons too.

What are the current practical limitations of explicitly correlated methods? What is the largest system that has been studied using explictly correlated methods?

Answer 1

Hylleraas wavefunctions tend only to be used for systems with up to 3 electrons. The reason for this is because the integrals have only been worked out analytically for up to 3 electrons, and they would be far too slow to do numerically. People have talked about using Hylleraas wavefunctions for 4 electrons, but when you ask them to show you results on real systems (for example the Be atom), the results are at least 12 orders of magnitude worse than what we have achieved for 3 electron systems, and about 4 orders of magnitude worse than what is achieved with methods based on non-Hylleraas wavefunctions.

Specifically, from Table I of this paper of mine on the carbon atom, you can see that the state-of-the art lowest variational energy for Be, is not obtained using a Hylleraas wavefunction but with an "explicitly correlated Gaussian". Nakatsuji has used Hylleraas-like wavefunctions for Be, but to the best of my knowledge, never achieved the correct energy to micro-Hartree precision, whereas explicitly correlated Gaussians have managed to achieve nano-Hartree precision almost a decade ago. Part of Table I from that paper is repeated here:

\begin{array}{c l l l} \textrm{He} & -2.903 724 377 034 119 598 311 159 245 194 404 446 696 925 309 838 & \textrm{Hylleraas-Log} & (2006)\\ \textrm{Li} & -7.478 060 323 910 134 843 & \textrm{Hylleraas} & (2017)\\ \textrm{Be} & -14.667 3564949 & \textrm{ECG} & (2013)\\ \end{array}

So while it is theoretically possible to use a 4-electron Hylleraas wavefunction, the necessary integrals have not yet been worked out analytically, and therefore no one has achieved any 4-electron energy to better accuracy than achieved using other methods (such as ECGs).

So now you may ask what the largest system is, that has been treated successfully using ECGs? ECGs have the advantage that the integrals can be calculated analytically. Unfortunately the answer is still 5 electrons, which was done by Puchalski et al. (the same author that did the Be atom in 2013) in 2015. I spoke to Puchalski many times about this, and while he did actually do some calculations on the 6-electron carbon atom, he did not obtain any results that he found to be worth publishing, although my 2018 paper on the carbon atom linked above (which used an aug-cc-pCV8Z basis set) inspired him to consider returning to the carbon atom project and publish something better than what I achieved (but that has not been completed yet). There has indeed been a 2019 paper which calculated the energy of the carbon atom using ECGs, but it did not achieve an energy lower than mine with the non-explicitly correlated aug-cc-pCV8Z basis set.

Nakatsuji has used Hylleraas-type wavefunctions for much larger numbers of electrons, but to the best of my knowledge, the results are not state-of-the-art for any of those systems.

Summary:

• Explicitly correlated Hylleraas wavefunctions have achieved state-of-the art energies up to 3e.
• Explicitly correlated Gaussian wavefunctions have achieved state-of-the art energies up to 5e.
  • However ECGs will likely soon achieve state-of-the-art energies for 6e$^-$.
  • ECGs will not likely achieve state-of-the-art energies for 7e$^-$ for several more years without a breakthrough change in the algorithm or hardware.
• Explcitly correlated wavefunctions do not tend to benefit from error cancellation as much as traditional single-particle basis sets like those from the Dunning family. Dunning prevails.

Answer 2

Hylleraas' method is special to the case of the helium atom; note that the equation only has a single nucleus. When you go to more electrons, you get more and more interelectronic distances in your wave function ansatz, which will blow up the scaling; this is why the method is impractical for several electrons. (Also, I think you would need several expansion parameters $\zeta$ to get an accurate result!)

The $r_{12}$ and $f_{12}$ methods are indeed inspired by Hylleraas' technique, but instead of aiming for the exact (i.e. FCI/CBS) energy their target is simpler: accelerate the basis set convergence for your approximate post-HF method like MP2. I would not say these methods are only rarely applied; you do see them quite often in e.g. benchmark studies. A rule of thumb you often hear is that these methods allow you to get one cardinal number more accurate results in your basis set, e.g. triple-$\zeta$ results for a double-$\zeta$ basis set, and quadruple-$\zeta$ results for a triple-$\zeta$ basis set. Since post-HF methods often have a steep scaling in the number of virtual orbitals (e.g. $O(v^4)$ for CCSD), this speedup can become important when you want to push the limits...