Could AAA-fit rationals replace Slater orbitals?

This page mostly exists because I forgot some basic algebra. I was reading about why gaussian basis functions were used and the main argument was that the 6d integral of two multi-center slater orbitals was not so fun. So I set out and brute-force integrated the functions. This was mostly proof that numerical integration of these things was a bad idea if you don't do it the smart way (NAOs). I was reading about pade approximants and realized — oh wow, rational functions can approximate a lot of functions so why not a slater-orbital. My ill-conceived idea was that you could turn the exponents into rational functions, multiply them together and then sum over the poles to do an integral. Very cool trick you can do with rational approximations. If you have a rational approximation you have an integral. I mean I guess it's pretty close to chebyshev interpolation in that respect, but I digress. The problem is that the centers are different still so the "multiply them together" trick doesn't work.

As an experiment, I ran a barycentric rational approximation using AAA, which can fit $e^{-r}$ to machine precision on a finite interval. So if nothing else you learn about AAA on this page.

AAA fits exp(−r) to machine precision on a bounded interval

Sample $e^{-r}$ on 4000 points across $r\in [0,20]$, run AAA with tolerance $10^{-13}$:

# AAA fit of exp(-r) on [0, 20], then test the most basic integral.
import numpy as np
from numpy.linalg import svd
from scipy.integrate import quad

def aaa(F, Z, tol=1e-13, mmax=100):
    """Greedy barycentric rational (Nakatsukasa-Sète-Trefethen 2018)."""
    Z, F = np.asarray(Z, complex), np.asarray(F, complex)
    M, J = len(Z), list(range(len(Z)))
    z, f = [], []
    C = np.zeros((M, 0), dtype=complex)
    R = np.mean(F) * np.ones(M, dtype=complex)
    for _ in range(mmax):
        jmax = J[int(np.argmax(np.abs(F[J] - R[J])))]
        z.append(Z[jmax]); f.append(F[jmax]); J.remove(jmax)
        with np.errstate(divide="ignore", invalid="ignore"):
            col = 1.0 / (Z - Z[jmax]); col[jmax] = 0.0
        C = np.column_stack([C, col])
        Cf = C[J, :]
        A = np.diag(F[J]) @ Cf - Cf @ np.diag(f)
        _, _, Vh = svd(A, full_matrices=False)
        w = Vh.conj().T[:, -1]
        num, den = C @ (w * np.array(f)), C @ w
        R = F.copy(); R[J] = num[J] / den[J]
        if J and np.max(np.abs(F[J] - R[J])) <= tol * np.max(np.abs(F)):
            break
    return np.array(z), np.array(f), np.array(w)

Z = np.linspace(0.0, 20.0, 4000)
z, f, w = aaa(np.exp(-Z), Z, tol=1e-13)

def Reval(r):
    t = w / (np.complex128(r) - z)
    return np.real(np.sum(t * f) / np.sum(t))

# Fit quality on the sample interval
rs = np.linspace(0.01, 20.0, 1001)
fit_err = max(abs(Reval(r) - np.exp(-r)) for r in rs)
print(f"AAA support points: {len(z)},   max fit error on [0, 20]: {fit_err:.2e}")

# Extrapolation beyond the sample interval
print("\nr        R(r)             exp(-r)")
for r in [21.0, 25.0, 30.0, 50.0, 100.0]:
    print(f"{r:>5.0f}   {Reval(r):+.3e}   {np.exp(-r):.3e}")

# The simplest STO integral: closed form is 1/4
I = quad(lambda r: r * r * Reval(r) ** 2, 0.0, np.inf, limit=400)[0]
print(f"\n∫₀ⁿ r² R(r)² dr  =  {I:.3e}    (closed form: 0.25)")

AAA support points: 11,   max fit error on [0, 20]: 4.95e-14

r        R(r)             exp(-r)
   21   +7.586e-10   7.583e-10
   25   +1.729e-10   1.389e-11
   30   +5.332e-09   9.358e-14
   50   +1.091e-06   1.929e-22
  100   +2.531e-05   3.720e-44

∫₀ⁿ r² R(r)² dr  =  2.401e+12    (closed form: 0.25)

Eleven support points get max fit error $5\times 10^{-14}$ on the whole sample interval. The convergence history is geometric: 10⁻¹ at iteration 2, 10⁻³ at iteration 4, 10⁻⁷ at iteration 8, machine precision at iteration 11. As a function approximation problem this is just done. The pole structure — 10 poles in complex-conjugate pairs at roughly 1 ± 15i, −3 ± 11i, −5 ± 8i, with appropriate residues — looks exactly like a high-order Padé approximant of an exponential, which is what AAA has effectively built.

Then the tail fails. Catastrophically.

The extrapolation table shows the failure directly. At $r=21$ the fit is still excellent — well outside the sample interval but still in the regime where polynomial behavior happens to track $e^{-r}$. A few units further out and the rational has no idea what an exponential is supposed to do. At $r=30$ it's already $6\times 10^4$ too big. At $r=100$ it's $10^{39}$ too big — the rational evaluates to $2.5\times 10^{-5}$, while $e^{-100}\approx 4\times 10^{-44}$.

This isn't a numerical glitch. It's structural. A rational function of degree $(m,n)$ with $m\le n$ decays at infinity like $r^{m-n}$ at best, i.e. polynomially. Exponential decay is unreachable by any finite rational. The Weierstrass approximation theorem only gives uniform approximation on compact sets; on an unbounded domain the rational can match the target arbitrarily well on the sample interval and remain completely uncontrolled outside it. AAA picked support points to drive down the fit error on $[0,20]$. It had no reason to behave at $r=50$, and it didn't.

The consequence: the most basic integral in quantum chemistry, $⟨\psi |\psi ⟩=4\pi {\int }_0^{\infty }{r}^2|\psi (r){|}^{2}dr$, blows up to $2\times 10^{12}$ instead of the closed-form $0.25$ (this is the radial part; the angular integration is trivial). The integrand $r^{2}R(r{)}^2$ is dominated entirely by the wrong-tail contribution from $r>30$, where AAA's rational is enormous and exp's is essentially zero. No amount of fit refinement on the sample interval fixes this — the problem is what the rational does where you weren't looking.

Classical [0/n] Padé saves integrability, loses fit quality

There's a different rational approximation that gives up AAA's accuracy but guarantees the right decay structure. The $[0/n]$ Padé approximant of $e^{-r}$ matches the Taylor series of $e^{-r}$ at $r=0$ through order $n$:

The numerator is constant (degree 0); the denominator has degree $n$. So $R_{[0/n]}(r)\sim n!/r^n$ at infinity by construction — guaranteed algebraic decay. The integral ${\int }_0^{\infty }{r}^2{R}^{2}dr$ converges as soon as $n\ge 2$.

# Classical [0/n] Padé of exp(-r):  R(r) = 1 / (1 + r + r²/2! + ... + rⁿ/n!).
# Numerator is constant => R decays like 1/rⁿ at infinity by construction.
import math, numpy as np
from scipy.integrate import quad

def pade_0n(n):
    Q = np.poly1d([1.0 / math.factorial(k) for k in range(n, -1, -1)])
    poles = np.array(Q.roots, dtype=complex)
    residues = np.array([1.0 / Q.deriv()(p) for p in poles])
    return poles, residues

def R_eval(r, poles, residues):
    return float(np.real(sum(c / (r - p) for c, p in zip(residues, poles))))

I_exact = 0.25
print(f"closed form  ∫ r² e⁻²ʳ dr  =  {I_exact}")
print(f"\n  n   max fit err on [0,10]      ∫ r² R² dr      Δ vs 0.25")
for n in [4, 6, 8, 10, 12]:
    poles, residues = pade_0n(n)
    rs = np.linspace(0.01, 10.0, 200)
    err = max(abs(R_eval(r, poles, residues) - np.exp(-r)) for r in rs)
    I = quad(lambda r: r * r * R_eval(r, poles, residues) ** 2, 0.0, np.inf, limit=400)[0]
    print(f" {n:>2}        {err:.2e}              {I:.10f}     {I - I_exact:+.2e}")

closed form  ∫ r² e⁻²ʳ dr  =  0.25

  n   max fit err on [0,10]      ∫ r² R² dr      Δ vs 0.25
  4        1.15e-02              0.2892171178    +3.92e-02
  6        2.29e-03              0.2551602527    +5.16e-03
  8        4.93e-04              0.2507580654    +7.58e-04
 10        1.10e-04              0.2501115925    +1.12e-04
 12        2.53e-05              0.2500161027    +1.61e-05

The integral converges geometrically to the closed-form $0.25$ as $n$ grows. But notice the fit-error column. At $n=12$, the max fit error on $[0,10]$ is $2.5\times 10^{-5}$. AAA hit $5\times 10^{-14}$ at $n=11$. Nine orders of magnitude difference — for fewer degrees of freedom on a similar interval. Classical Padé is locking down the Taylor behavior at $r=0$ instead of distributing error optimally across the interval, which AAA does freely. The cost of guaranteeing the tail decay is throwing away most of the fit quality.

The one positive result: residue summation does work, on the integrable rational

With the $[0/n]$ Padé in pole-residue form $R(r)={\sum }_k{c}_k/(r-p_k)$ with ${\sum }_k{c}_k=0$, we can compute ${\int }_0^{\infty }{r}^{2}R(r{)}^{2}dr$ analytically via partial-fraction decomposition of $R(r{)}^2$:

# Same [0/n] Padé, but the integral done analytically via residue summation
# instead of scipy.quad. Building blocks:
#
#   J(L; p, q) = ∫₀ᴸ r² / [(r-p)(r-q)] dr  for p ≠ q
#              = L + ( p² log((L-p)/(-p)) - q² log((L-q)/(-q)) ) / (p - q)
#
#   K(L; p)    = ∫₀ᴸ r² / (r-p)² dr
#              = L + 2p log((L-p)/(-p)) + p² (-1/p - 1/(L - p))
#
# Each (k != j) cross term contributes c_k c_j J(L; p_k, p_j).
# Each (k = j) diagonal term contributes c_k² K(L; p_k).
# Leading L pieces cancel across the double sum because Σ c_k = 0.
import math, numpy as np

def J_offdiag(L, p, q):
    return L + (p * p * np.log((L - p) / (-p)) - q * q * np.log((L - q) / (-q))) / (p - q)

def K_diag(L, p):
    return L + 2 * p * np.log((L - p) / (-p)) + p * p * (-1.0 / p - 1.0 / (L - p))

def integral_via_residues(L, poles, residues):
    total = 0.0 + 0.0j
    for k, (ck, pk) in enumerate(zip(residues, poles)):
        for j, (cj, pj) in enumerate(zip(residues, poles)):
            c = ck * cj
            total += c * (K_diag(L, pk) if k == j else J_offdiag(L, pk, pj))
    return total.real

# pade_0n() as above
print("  n   residue sum, L=500    residue sum, L=2000   Δ vs scipy.quad")
for n in [4, 6, 8, 10, 12]:
    poles, residues = pade_0n(n)
    I_500 = integral_via_residues(500.0, poles, residues)
    I_2k  = integral_via_residues(2000.0, poles, residues)
    print(f" {n:>2}      {I_500:.10f}         {I_2k:.10f}      {I_2k - I_quad:+.2e}")

  n   residue sum, L=500    residue sum, L=2000   Δ vs scipy.quad
  4      0.2892171178         0.2892171178      -2.71e-12
  6      0.2551602527         0.2551602527      +8.13e-12
  8      0.2507580654         0.2507580654      +8.32e-12
 10      0.2501115925         0.2501115926      +1.20e-10
 12      0.2500161027         0.2500161028      +3.57e-11

The leading-$L$ pieces cancel across the double sum because ${\sum }_k{c}_k=0$ (forced by the decay structure of a strictly proper rational). What's left is a finite expression in pole locations and residues that agrees with scipy.quad to $10^{-11}$. The cutoff $L$ matters in principle — log terms have an $L$-dependent piece — but their contribution drops as $1/L$ and the result is essentially exact at $L=500$ already. So the pole-summation idea works once the rational has the right decay structure. The piece of the original question that survives.

The second wall: multi-center, and where the rational idea actually breaks

Suppose we fix wall one (use $[0/n]$ Padé instead of AAA, accept the loss of fit accuracy) and put two such basis functions on different centers $\vec{A}$ and $\vec{B}$. The overlap integral $\int {\chi }_A(|\vec{r}-\vec{A}|){\chi }_B(|\vec{r}-\vec{B}|)d^{3}r$ has integrand $R_A(|\vec{r}-\vec{A}|)\cdot R_B(|\vec{r}-\vec{B}|)$ — a product of two rationals in two different distances from two different centers.

Here's where the original motivation goes wrong, and it's worth naming clearly. The reason I thought rationals might help in the first place was: with Slater orbitals, the exponential $e^{-\zeta |\vec{r}-\vec{A}|}$ traps the distance $|\vec{r}-\vec{A}|$ inside a transcendental function. Multiplying two such exponentials on different centers gives $\exp (-{\zeta }_A|\vec{r}-\vec{A}|-{\zeta }_B|\vec{r}-\vec{B}|)$, and that exponent doesn't factor — the sum of two distances has no nice closed form. Pull the distance out of the exponential by replacing it with a rational, and now $|\vec{r}-\vec{A}|$ is algebraically accessible — you can take its poles, do partial fractions, manipulate it.

But the difficulty was never that the variable was stuck inside an exponential. The difficulty is that $|\vec{r}-\vec{A}|$ and $|\vec{r}-\vec{B}|$ are two different functions of the integration variable $\vec{r}$. Pulling them out of the exponential doesn't make them the same function. The product $R_A(|\vec{r}-\vec{A}|)\cdot R_B(|\vec{r}-\vec{B}|)$ still has the two-center geometry baked in: partial fractions in $|\vec{r}-\vec{A}|$ leaves $|\vec{r}-\vec{B}|$ untouched, and vice versa. The integral doesn't reduce.

What's actually special about Gaussians is that they're not just exponentials — they're exponentials of quadratic functions of $\vec{r}$. The sum $\alpha |\vec{r}-\vec{A}{|}^2+\beta |\vec{r}-\vec{B}{|}^2$ is quadratic in $\vec{r}$, and any quadratic completes the square:

with $\vec{P}=(\alpha \vec{A}+\beta \vec{B})/(\alpha +\beta )$. A sum of two quadratics in $\vec{r}$ is a single quadratic in $\vec{r}$ plus a constant. That's what factors. Exponentiate both sides and you get the Boys product theorem — two Gaussians on different centers fuse into one Gaussian on a third center, times a closed-form prefactor. The exponential isn't doing the work; the Pythagorean identity on the quadratic underneath it is.

Slater orbitals depend on distance linearly inside the exponent, so they have no completing-the-square identity. Rationals depend on distance through a polynomial-ratio, so they have none either. Whatever you replace the Gaussian with, the algebraic identity that makes 4-center 2-electron integrals tractable is specific to quadratic dependence on $\vec{r}$. The rational idea wasn't a different approach to the same problem; it was a different approach to the wrong problem. Solving "the exponential is the obstruction" doesn't help if the actual obstruction is "two different distances don't simplify together."

The actual lesson: rationals work inside exponentials, never instead of them

The original motivation was that rational functions have a tractable pole-residue structure. They do. But the way that structure succeeds in electronic structure is by sitting inside an exponential, not by replacing one. The Padé-Jastrow factor is the canonical example:

A Slater determinant $D$ multiplied by an exponential of a sum of rational two-body terms. The rational shape captures the electron-electron cusp (it's linear in $r$ near zero, with slope fixed by the Kato cusp condition: ${\omega }_0=1/2$ for antiparallel spins, $1/4$ for parallel). The exponential wrapping makes the whole wavefunction normalizable and integrable over ${\mathbb{R}}^{3N}$. Neither factor works on its own; each fixes what the other can't do.

The $\exp$ isn't there for looks: it fixes the tail (exponential decay of the full wavefunction), it makes the Jastrow factor strictly positive (so the wavefunction's nodal structure is controlled by $D$ alone), and it puts the whole thing in a form that QMC can sample efficiently ($|\Psi {|}^2=e^{2J}|D{|}^2$ stays positive, so importance sampling works). Strip the exp out and use the rational as $\psi$ directly — like the variational rational page does — and you lose all three properties.

Where pole expansions DO win in electronic structure

The pole-residue trick has real wins in the field, just not for basis functions. The common thread: it wins where the function being approximated is intrinsically rational (or has a natural contour representation), not where it's being forced into a rational form against its nature.

• PEXSI (pole expansion + selected inversion). Lin Lin et al. approximate the Fermi-Dirac function by a sum of complex poles (typically 30-80), then evaluate density-matrix elements via selected inversion of sparse matrices at those poles. Replaces $O(N^3)$ diagonalization with $O(N)$ sparse work for the right class of systems. The Fermi-Dirac function is intrinsically meromorphic, so a finite rational fit is the natural representation.
• GW with contour deformation. The self-energy involves a frequency integral that has slowly-decaying tails on the real axis but well-behaved poles in the complex plane. Deform the contour, sum residues. Same idea, different operator.
• AAA-fit resolvents for Krylov methods. To compute $(zI-A{)}^{-1}b$ at many shifts $z$, AAA-fit the resolvent matrix element $b^T(zI-A{)}^{-1}b$ in $z$ from samples produced by a Lanczos chain, then evaluate the rational at any shift cheaply. Here the function in $z$ is genuinely rational — it's the Padé-via-Lanczos continued fraction — so AAA is recovering a structural property rather than imposing one.

In every one of these cases the rational form is matching the underlying mathematics. Basis functions in real space aren't rational — they want to be exponentials with a cusp — so dressing them in a rational doesn't recover any structure that was there to begin with. The exercise on this page makes the failure mode concrete: the rational can be made arbitrarily accurate on a bounded sample interval, or made to have integrable tails, but not both.

What I'd still try, given infinite spare time

Three half-formed follow-ups worth listing, in case any of them actually works:

1. Variable transform. Map $r\in [0,\infty )$ to $t\in [0,1)$ via $t=r/(1+r)$, fit $e^{-r(t)}$ by AAA in $t$. The fit is now on a compact interval so the tail problem goes away. But the basis function is no longer rational in $r$, only in $t$ — and the volume element $r^{2}dr$ becomes $t^2/(1-t{)}^{4}dt$, which breaks the partial-fraction structure when integrated against another $t$-rational on a different center. Probably another negative result, but worth checking cleanly.
2. Exp(rational) directly, not Jastrow-style. Use $\psi (r)=\exp (-\alpha r+g(r))$ where $g$ is a small rational correction. Now you have an explicit single-particle Slater-Jastrow form. The radial integrals can be done because the exponential decay is preserved by the linear-in-r leading behavior. The 4-center 2-electron integral still doesn't reduce analytically, but for small atoms a sparse-pole correction to Gaussians might give a useful pre-Jastrow improvement at low cost. This is approximately what Maximally-Compact Gaussian-based methods already do in spirit.
3. Use the AAA fit only where it's trustworthy. If AAA gives you a rational that's accurate on $[0,R_{\max }]$ and garbage outside it, restrict the integration domain to $[0,R_{\max }]$ with a windowing function. Combine with a decaying tail correction (e.g. continue the integration with an asymptotic exp form). This is hacky but might be a practical way to use AAA for STO-like profiles in a controlled regime — relevant for basis-set construction more than for direct integration.

Related on this site

• Variational rational trial on hydrogen — the pedagogical case of using a rational as $\psi$ directly. Shows the bound degrades but stays finite, because there the integration was set up correctly.
• Variational principle — exponential trial baseline, which gets $E_0=-0.5$ on the nose.
• Padé approximants — what [m/n] Padé means structurally, and where it works for resumming divergent series.
• Lanczos iteration — the Padé-via-Lanczos correspondence that makes AAA-fit resolvents work where AAA-fit basis functions don't.