Basis Sets: Slater vs Gaussian Functions

Every quantum chemistry calculation begins with a lie. That lie turns out to be extraordinarily useful. Instead of solving for the true, unbounded shape of a molecular orbital, we approximate it. We use a sum of simple, pre-chosen mathematical functions. The choice of those functions — the basis set — determines whether a calculation runs in seconds. Or whether it never finishes at all. This article is about two candidate shapes. They have shaped the entire history of computational chemistry. The first is the Slater-type orbital (STO). It mimics the true hydrogen-atom wavefunction almost perfectly. The second is the Gaussian-type orbital (GTO). It is a physically worse approximation but an algebraically miraculous one. The stakes are not academic. The Roothaan-Hall Equations reduce the Hartree-Fock problem to a matrix eigenvalue problem. But this only works once every orbital is expanded in some finite basis. That expansion requires computing enormous numbers of two-electron repulsion integrals. For a system with $N$ basis functions, there are of order $N^4$ such integrals. If each one requires a slow numerical quadrature, even modest molecules become computationally hopeless. The entire practical existence of ab initio quantum chemistry rests on one fact. This includes everything from a Hartree-Fock calculation on water to a density functional theory geometry optimization on a drug candidate. The fact is this: a specific class of integrals over Gaussian functions can be evaluated in closed form. This can be done quickly and exactly. The big picture is a trade. Slater functions get the physics right. They capture the correct cusp at the nucleus. They also capture the correct exponential decay at long range. Both are dictated by the shape of true atomic orbitals. But Slater functions make the integrals nearly intractable for polyatomic molecules. Gaussian functions get the physics slightly wrong at both extremes. However, the product of two Gaussians is another Gaussian. This fact is known as the Gaussian Product Theorem. This single algebraic fact converts multicenter integrals into single-center ones. Those single-center integrals can be solved analytically. Modern software resolves the tension by fitting linear combinations of Gaussians. These combinations reproduce Slater-like shapes. This buys back most of the missing physics. It also keeps all of the computational speed. By the end of this article, you will be able to write down both orbital forms. You will derive the Gaussian Product Theorem from first principles. You will understand why STO-based integrals defeated early quantum chemists. You will also see exactly how contracted Gaussian basis sets solve the problem in practice. This sets the stage for the natural next question. How many Gaussians, and of what type, do you actually need per atom? That is the subject of minimal and split-valence basis set design.

Key takeaways

Slater-type orbitals (STOs) have the form $e^{-\zeta r}$ , matching the true nuclear cusp and long-range decay of atomic orbitals exactly.
Gaussian-type orbitals (GTOs) have the form $e^{-\alpha r^2}$ : wrong at $r=0$ (no cusp) and wrong at large $r$ (decays too fast), but algebraically tractable.
The Gaussian Product Theorem states that a product of two Gaussians on different centers is a single Gaussian on a third center — this is the reason GTOs win.
Multicenter two-electron integrals over GTOs reduce to closed-form expressions (the Boys function); the same integrals over STOs generally require numerical multicenter quadrature.
Practical basis sets use contracted Gaussians — fixed linear combinations of several primitive GTOs — to approximate the STO shape while keeping GTO speed.
Nearly all modern electronic structure codes (Gaussian, ORCA, Psi4, NWChem, Q-Chem) use contracted GTO basis sets for exactly this reason.
Key idea: GTOs are a worse model of an atomic orbital than STOs, but their algebraic closure under multiplication makes molecular integrals analytically solvable, and that computational advantage decides which function wins in practice.

Where this lives today

Every mainstream quantum chemistry package — Gaussian, ORCA, Psi4, NWChem, Q-Chem, Molpro, and the Gaussian-basis modules of many periodic DFT codes — represents molecular orbitals as linear combinations of contracted Gaussian-type basis functions for exactly the reason developed here: the Roothaan-Hall Equations require the Fock matrix's two-electron repulsion integrals to be computed millions to billions of times in a single geometry optimization or ab initio molecular dynamics run, and only Gaussian integrals stay analytically tractable at that scale. Drug-discovery workflows that use DFT to refine docking poses, materials-science calculations of band gaps in candidate photovoltaic crystals, and high-accuracy coupled-cluster benchmarks of reaction barrier heights in atmospheric chemistry all inherit this same computational bedrock, whether or not the practitioner ever thinks about Boys functions directly. Slater-type orbitals have not disappeared; they persist where their superior physical shape outweighs the integral cost, most notably in the Amsterdam Density Functional (ADF) program, which uses Slater-type basis functions for molecular and periodic DFT by exploiting numerical integration grids rather than fully analytical multicenter integrals, and in semi-empirical methods (extended Hückel theory, early PPP and CNDO/INDO schemes) where STOs' compactness and correct decay were prized and expensive multicenter integrals were approximated or discarded rather than computed exactly. The contraction strategy foreshadowed in Example 2 — fitting several Gaussian primitives to approximate an STO shape, then freezing the coefficients — is the direct ancestor of every basis set name a computational chemist encounters daily: STO-3G, 6-31G, cc-pVDZ, def2-TZVP. Understanding why those names contain the word "Gaussian" implicitly, and why their designers still chase STO-like shapes explicitly, is the natural next step, taken up in the study of minimal and split-valence basis set construction.

Ask the tutor: How do STO-3G, 6-31G, and cc-pVDZ differ in how many contracted Gaussians they assign to core versus valence orbitals?

The physical picture

An atomic orbital in a many-electron atom is not exactly hydrogenic, but it retains two structural features inherited from the exact hydrogen-atom solutions. Near the nucleus, the electron feels a $-Z/r$ Coulomb singularity, and the Schrödinger equation forces the wavefunction to develop a cusp there — a kink with slope proportional to $-Z$ at $r=0$ (Kato's cusp condition). Far from the nucleus, in the classically forbidden region, the wavefunction decays as a pure exponential, $e^{-\zeta r}$ , with $\zeta$ related to the ionization energy of that electron. Both features are baked into the Slater-type orbital by construction:

\chi^{\text{STO}}_{n\ell m}(r,\theta,\phi) = N \, r^{n-1} e^{-\zeta r} \, Y_{\ell m}(\theta,\phi)

This is essentially a hydrogen-like orbital with the principal quantum number's polynomial prefactor simplified and the effective nuclear charge folded into a single adjustable parameter $\zeta$ (the Slater exponent). STOs are, by design, an excellent one-center representation of an atomic orbital: correct cusp, correct tail, few nodes required to get the shape right. The Gaussian-type orbital instead uses $e^{-\alpha r^2}$ as the radial decay:

\chi^{\text{GTO}}_{\ell m}(x,y,z) = N \, x^{l} y^{m} z^{n} e^{-\alpha r^2}

(written here in the common Cartesian form, with $l+m+n$ giving the total angular momentum). A single GTO is a poor model of an orbital in two specific, well-understood ways. At $r=0$ , $d(e^{-\alpha r^2})/dr = 0$ — the function is smooth at the nucleus, with zero slope, whereas the true wavefunction has a sharp cusp. At large $r$ , $e^{-\alpha r^2}$ falls off much faster than $e^{-\zeta r}$ , so a single GTO underestimates the long-range tail of the electron density that governs, for instance, van der Waals interactions and diffuse anion orbitals. Neither defect is small: plotted together, a single STO and a single GTO with exponents tuned to match at one point diverge badly everywhere else. Yet the entire field of computational chemistry chose GTOs, and the reason is purely algebraic, not physical: the Gaussian Product Theorem. The product of two 1-dimensional (or 3-dimensional) Gaussians centered at two different nuclei $A$ and $B$ is itself a single Gaussian, centered at a point on the line between them, with a combined exponent and a multiplicative prefactor that depends only on the separation $|A-B|$ . No such closure property holds for STOs: the product of two decaying exponentials centered at different points is not a simple exponential centered anywhere convenient, and multicenter integrals involving STOs generically require elliptic-coordinate expansions or numerical quadrature that scale poorly as the number of atoms grows. Because Hartree-Fock and post-Hartree-Fock methods require the two-electron repulsion integral

(\mu\nu|\lambda\sigma) = \iint \chi_\mu(1)\chi_\nu(1)\frac{1}{r_{12}}\chi_\lambda(2)\chi_\sigma(2)\, d\tau_1 d\tau_2

over basis functions that in general sit on four different atomic centers, the Gaussian Product Theorem is not a mathematical curiosity — it is the mechanism that makes molecular quantum chemistry computationally feasible at all. Applying the theorem twice (once to $\chi_\mu\chi_\nu$ , once to $\chi_\lambda\chi_\sigma$ ) collapses the four-center integral into a two-center Coulomb integral between two Gaussian charge distributions, which itself reduces to a well-tabulated special function (the Boys function) after one more application of Gaussian integral identities.

Figure

Cusp and tail: STO vs GTO radial shapes near and far from the nucleus

Figure

Gaussian Product Theorem: two Gaussians on centers A, B collapse to one Gaussian on center P

The core ideas

The formal statement of the Gaussian Product Theorem, restricted for clarity to one Cartesian dimension (the three-dimensional case factors identically in $x$ , $y$ , $z$ ), is:

e^{-\alpha (x-A)^2} \, e^{-\beta (x-B)^2} = K_{AB}\; e^{-\gamma (x-P)^2}

where:

$\alpha, \beta$ are the exponents of the two original Gaussians, centered at $A$ and $B$ ,
$\gamma = \alpha + \beta$ is the exponent of the resulting single Gaussian,
$P = \dfrac{\alpha A + \beta B}{\alpha + \beta}$ is the new center, a weighted average of $A$ and $B$ ,
$K_{AB} = \exp\!\left[-\dfrac{\alpha\beta}{\alpha+\beta}(A-B)^2\right]$ is a constant prefactor that depends only on the separation $A-B$ , not on $x$ .

Reading: no matter how far apart $A$ and $B$ are, multiplying two Gaussians always yields exactly one Gaussian at a new, easily computed center, times a number you can compute immediately. This is not an approximation — it is an algebraic identity valid for every $\alpha, \beta, A, B$ . The practical payoff appears in the two-electron repulsion integral. Because a product of basis functions on two different centers becomes one Gaussian "charge overlap distribution," the four-center integral

(\mu\nu|\lambda\sigma)

reduces, after applying the theorem to the $\mu\nu$ pair and the $\lambda\sigma$ pair, to an integral between two single Gaussian charge clouds centered at $P$ and $Q$ . That remaining two-center Coulomb integral has a closed analytical form involving the Boys function:

F_0(T) = \int_0^1 e^{-T t^2}\, dt = \frac{1}{2}\sqrt{\frac{\pi}{T}}\,\mathrm{erf}\!\left(\sqrt{T}\right), \qquad T = \gamma_{PQ}\,|P-Q|^2

where:

$T$ is a dimensionless argument built from the combined exponent of the $P$ -Gaussian and $Q$ -Gaussian pair and the squared distance between their centers,
$\mathrm{erf}$ is the standard error function, itself expressible via elementary functions or rapid series/asymptotic expansions,
higher angular momentum integrals require $F_n(T)$ for $n>0$ , obtained from $F_0$ by a stable downward or upward recursion.

Reading: every single electron repulsion integral over $s$ -type Gaussians reduces to evaluating one Boys function at one number $T$ . This is the computational miracle: a 4-dimensional, 4-center, singular integral has become the evaluation of a tabulated special function. No comparably simple closed form exists for the analogous integral over four Slater-type orbitals sitting on four different nuclei; those integrals in general require numerical evaluation in prolate spheroidal (elliptic) coordinates, feasible for diatomics but prohibitively slow for the multicenter case of an ordinary molecule.

Derivation AThe Gaussian Product Theorem, from a completed square

Goal. Show algebraically that $e^{-\alpha(x-A)^2}e^{-\beta(x-B)^2}$ is a single Gaussian in $x$ , and identify its center, exponent, and prefactor.

Step 1 — Expand both exponents. Write out the two quadratics explicitly:

-\alpha(x-A)^2 - \beta(x-B)^2 = -\alpha x^2 + 2\alpha A x - \alpha A^2 - \beta x^2 + 2\beta B x - \beta B^2

Collect the terms by power of $x$ :

= -(\alpha+\beta)x^2 + 2(\alpha A + \beta B)x - (\alpha A^2 + \beta B^2)

Step 2 — Define the combined exponent and center. Let $\gamma \equiv \alpha + \beta$ and define

P \equiv \frac{\alpha A + \beta B}{\alpha + \beta}

so that the linear-in- $x$ coefficient is $2\gamma P$ . The exponent becomes

-\gamma x^2 + 2\gamma P x - (\alpha A^2 + \beta B^2)

Step 3 — Complete the square in $x$ . Add and subtract $\gamma P^2$ :

-\gamma x^2 + 2\gamma P x - \gamma P^2 + \gamma P^2 - (\alpha A^2 + \beta B^2) = -\gamma(x-P)^2 + \left[\gamma P^2 - \alpha A^2 - \beta B^2\right]

The first term is exactly the exponent of a single Gaussian centered at $P$ with exponent $\gamma$ .

Step 4 — Simplify the leftover constant. Substitute $\gamma P = \alpha A + \beta B$ and $\gamma = \alpha+\beta$ into the bracket:

\gamma P^2 - \alpha A^2 - \beta B^2 = P(\alpha A + \beta B) - \alpha A^2 - \beta B^2 = \alpha A(P-A) + \beta B(P - B)

Now use $P - A = \dfrac{\beta(B-A)}{\alpha+\beta}$ and $P - B = \dfrac{\alpha(A-B)}{\alpha+\beta}$ (direct from the definition of $P$ ):

= \frac{\alpha\beta A(B-A)}{\alpha+\beta} + \frac{\alpha\beta B(A-B)}{\alpha+\beta} = \frac{\alpha\beta(A-B)}{\alpha+\beta}\left[B - A\right]\cdot(-1)\cdots = -\frac{\alpha\beta}{\alpha+\beta}(A-B)^2

(the algebra collapses cleanly because the two terms are proportional to $(A-B)$ with opposite signs on $A$ and $B$ , giving a single $-(A-B)^2$ term).

Step 5 — Assemble the result.

e^{-\alpha(x-A)^2}e^{-\beta(x-B)^2} = \exp\left[-\frac{\alpha\beta}{\alpha+\beta}(A-B)^2\right]\exp\left[-\gamma(x-P)^2\right]

which is precisely the theorem stated above, with $K_{AB} = \exp\left[-\frac{\alpha\beta}{\gamma}(A-B)^2\right]$ . Each of the three Cartesian directions factors identically, so the full 3-D theorem follows by taking the product over $x,y,z$ .

Go deeper · Derivation B · 2 minWhy the same closure fails for STOs, and what it costs

Derivation BWhy the same closure fails for STOs, and what it costs

Goal. Show that the product of two Slater-type exponentials on different centers is not reducible to a single one-center exponential, and identify the computational consequence.

Step 1 — Write the STO product in the naive way. Consider two 1s-like STOs on centers $A$ and $B$ :

e^{-\zeta_1 r_A}\, e^{-\zeta_2 r_B}, \qquad r_A = |\mathbf r - \mathbf A|,\ r_B = |\mathbf r - \mathbf B|

Because $r_A$ and $r_B$ are square roots of quadratics in $\mathbf r$ (not quadratics themselves), there is no algebraic completion-of-the-square step analogous to Step 3 above. The sum $\zeta_1 r_A + \zeta_2 r_B$ does not collapse into $\zeta_{\text{eff}}\, |\mathbf r - \mathbf P|$ for any fixed point $P$ and constant $\zeta_{\text{eff}}$ , because the level surfaces of $r_A + r_B$ (constant) are prolate spheroids (ellipsoids of revolution about the $A$ – $B$ axis), not spheres about any single center.

Step 2 — Consequence for two-center integrals. For a diatomic molecule this is still tractable: switching to prolate spheroidal coordinates $\mu = (r_A+r_B)/R$ , $\nu=(r_A-r_B)/R$ (with $R=|A-B|$ ) diagonalizes the geometry, and overlap, kinetic, and even some two-electron integrals over STOs can be done via tabulated auxiliary functions ( $A_n(\mu)$ , $B_n(\nu)$ integrals). This is exactly how early diatomic STO Hartree-Fock calculations (e.g., on $\mathrm{N_2}$ , $\mathrm{CO}$ ) were performed by hand and by early computers.

Step 3 — Consequence for multicenter integrals. A general two-electron integral $(\mu\nu|\lambda\sigma)$ in a polyatomic molecule involves basis functions on up to four different nuclei simultaneously. There is no coordinate system in which all four Slater exponentials simultaneously become separable, single-center functions. The integral must be attacked by multipole/Fourier-transform expansions or brute numerical quadrature in 6 dimensions (3 for each electron), and convergence is slow because of the STO's nucleus-centered cusp, which itself is hard to represent on a numerical grid.

Step 4 — Scaling consequence. Because Gaussian integrals reduce to closed forms independent of the number of centers involved, GTO-based codes can evaluate $(\mu\nu|\lambda\sigma)$ integrals with the same low computational cost (dominated by evaluating a Boys function and simple prefactors) regardless of whether $\mu,\nu,\lambda,\sigma$ sit on one, two, three, or four different atoms. STO-based codes pay an escalating cost as the number of distinct centers per integral increases — the very situation that occurs constantly in polyatomic molecules with dozens or hundreds of atoms.

Step 5 — The engineering resolution. Rather than abandon the STO's superior physical shape, chemists fit a fixed linear combination of several Gaussian primitives to reproduce an STO's radial profile as closely as possible:

\chi^{\text{STO-}n\text{G}}(r) = \sum_{k=1}^{n} d_k\, \left(\frac{2\alpha_k}{\pi}\right)^{3/4} e^{-\alpha_k r^2}

where the coefficients $d_k$ and exponents $\alpha_k$ are pre-optimized (e.g., by least-squares fit to the STO or by direct energy minimization) and then frozen — this is a contracted Gaussian basis function. Every integral is still computed over primitive Gaussians (retaining full analytic tractability), but the resulting contracted function behaves far more like a true STO, especially near the nucleus and in its overall radial envelope, than any single Gaussian could.

Worked examples

Example 1 — Foundational: applying the Gaussian Product Theorem to two off-center $s$ -Gaussians. Take two normalized primitive $s$ -type Gaussians with unequal exponents, placed on two different nuclei separated along the $x$ -axis: $\alpha = 0.50\ \text{bohr}^{-2}$ centered at $A = 0.00$ bohr, and $\beta = 1.20\ \text{bohr}^{-2}$ centered at $B = 1.40$ bohr. (Deliberately asymmetric numbers, so the formula is genuinely exercised rather than trivially satisfied.) The combined exponent is

\gamma = \alpha + \beta = 0.50 + 1.20 = 1.70\ \text{bohr}^{-2}

The new center is the exponent-weighted average:

P = \frac{\alpha A + \beta B}{\alpha+\beta} = \frac{(0.50)(0.00) + (1.20)(1.40)}{1.70} = \frac{1.68}{1.70} = 0.9882\ \text{bohr}

Note $P$ sits much closer to $B$ than to $A$ , exactly as expected since $\beta > \alpha$ — the larger exponent (more sharply localized function) pulls the product center toward itself. The prefactor requires $(A-B)^2 = (0.00 - 1.40)^2 = 1.96\ \text{bohr}^2$ :

K_{AB} = \exp\left[-\frac{\alpha\beta}{\alpha+\beta}(A-B)^2\right] = \exp\left[-\frac{(0.50)(1.20)}{1.70}(1.96)\right]

= \exp\left[-\frac{0.60}{1.70}\times 1.96\right] = \exp\left[-0.3529\times 1.96\right] = \exp[-0.6918] = 0.5007

So the product collapses exactly to

e^{-0.50(x-0)^2}\,e^{-1.20(x-1.40)^2} = 0.5007\; e^{-1.70(x-0.9882)^2}

Limiting-case check: if the two centers coincide, $A=B$ , then $(A-B)^2=0$ so $K_{AB}=1$ and $P=A=B$ trivially — the product of two Gaussians on the same center is just a single Gaussian of exponent $\gamma=\alpha+\beta$ on that center, with no attenuation, which is obviously correct by direct multiplication of exponents. As the separation $|A-B|\to\infty$ with fixed $\alpha,\beta$ , $K_{AB}\to 0$ exponentially — the overlap of two well-separated Gaussians vanishes, again matching physical intuition (electrons on distant atoms barely overlap).

Example 2 — Applied: comparing STO and single-GTO representations of hydrogen 1s, and the "best-fit" GTO exponent. The exact hydrogen 1s orbital is an STO with $\zeta = 1.00\ \text{bohr}^{-1}$ :

\chi^{\text{STO}}_{1s}(r) = \left(\frac{1}{\pi}\right)^{1/2} e^{-r}

If we instead use a single normalized $s$ -Gaussian $\chi^{\text{GTO}}(r) = (2\alpha/\pi)^{3/4} e^{-\alpha r^2}$ , the exponent $\alpha$ that minimizes the Hartree-Fock ground-state energy of the hydrogen atom (equivalently maximizes overlap while respecting the variational principle) is the well-known value

\alpha^{*} = \frac{8}{9\pi} \approx 0.2829\ \text{bohr}^{-2}

Using this single "best" Gaussian, the resulting hydrogen-atom energy is

E_{\text{1-GTO}} = -\frac{4}{3\pi} \; \text{hartree} \approx -0.4244\ \text{hartree}

Compare this with the exact energy $E_{\text{exact}} = -0.5000$ hartree. The single Gaussian is in error by

\Delta E = -0.4244 - (-0.5000) = +0.0756\ \text{hartree} \approx +47.5\ \text{kcal/mol}

This error — over one Hartree unit's worth of chemical accuracy lost from a single misplaced function shape — is the quantitative face of the missing cusp: near $r=0$ the Gaussian's flat top places electron density in the wrong place relative to the sharply peaked exact function, and the variational energy pays the price. Now contract several Gaussians together (an "STO- $n$ G" fit, in the spirit of the derivation above). Using the standard STO-3G fit to $\zeta=1.0$ (exponents $\alpha_1=0.1688554$ , $\alpha_2=0.6239137$ , $\alpha_3=3.425251$ bohr $^{-2}$ , appropriately weighted and normalized), the resulting ground-state energy recovers to

E_{\text{STO-3G}} \approx -0.4666\ \text{hartree}

closing more than half the gap to the exact $-0.5000$ hartree using just three Gaussian primitives contracted into one function — while every integral in the calculation remains a fast, closed-form Gaussian integral. Adding more primitives (STO-6G, and beyond) closes the remaining gap essentially completely, at the cost of more primitive integrals per contracted function, not a change in integral type.

Example 3 — Exam/competition level: the Boys function argument for a heteronuclear overlap-derived Coulomb integral. Consider two normalized $s$ -Gaussian charge distributions formed from Gaussian products (as in Example 1): one centered at $P_1 = 0.30$ bohr with combined exponent $\gamma_1 = 0.90\ \text{bohr}^{-2}$ , the other centered at $P_2 = 2.10$ bohr with combined exponent $\gamma_2 = 1.50\ \text{bohr}^{-2}$ . We want the argument $T$ of the Boys function that enters their Coulomb repulsion integral (the standard combination that appears once the two-electron integral has been reduced by the product theorem on each electron separately):

T = \frac{\gamma_1 \gamma_2}{\gamma_1+\gamma_2}\,|P_1 - P_2|^2

First the reduced exponent:

\frac{\gamma_1\gamma_2}{\gamma_1+\gamma_2} = \frac{(0.90)(1.50)}{0.90+1.50} = \frac{1.35}{2.40} = 0.5625\ \text{bohr}^{-2}

Then the squared separation:

|P_1-P_2|^2 = (2.10-0.30)^2 = (1.80)^2 = 3.24\ \text{bohr}^2

T = 0.5625 \times 3.24 = 1.8225

The Boys function at this argument,

F_0(1.8225) = \frac{1}{2}\sqrt{\frac{\pi}{1.8225}}\,\mathrm{erf}\!\left(\sqrt{1.8225}\right)

requires $\sqrt{1.8225} = 1.3500$ and $\mathrm{erf}(1.3500) \approx 0.9438$ (from standard error-function tables). Then

\sqrt{\pi/1.8225} = \sqrt{1.7237} = 1.3129

F_0(1.8225) = \frac{1}{2}(1.3129)(0.9438) = 0.6197

Limiting case check: as $T\to 0$ (the two charge clouds collapse onto the same point, or the exponents both shrink toward zero making both distributions infinitely diffuse and overlapping), $F_0(T) \to 1$ , since $\mathrm{erf}(\sqrt T)/\sqrt T \to 2/\sqrt\pi$ in that limit and the prefactors cancel exactly — this is the correct "maximum overlap, minimum separation" limit, matching the known series expansion $F_0(T) = 1 - T/3 + T^2/10 - \cdots$ for small $T$ . As $T\to\infty$ (well-separated, sharply localized charges), $F_0(T)\to \frac{1}{2}\sqrt{\pi/T}\to 0$ , consistent with vanishing wavefunction overlap at large separation. Both checks confirm the closed form behaves exactly as electrostatics demands, and this single number $F_0(1.8225)=0.6197$ — obtained without any numerical integration — is precisely what a Gaussian-based Hartree-Fock program computes millions of times per SCF cycle on a real molecule.

Limiting cases, scope & common traps

Regime / question	STO behavior	GTO behavior
At the nucleus ( $r\to0$ )	Correct cusp, $d\chi/dr\|_{0}\ne 0$	No cusp, $d\chi/dr\|_0 = 0$ (must use many primitives to approximate)
Long range ( $r\to\infty$ )	Correct $e^{-\zeta r}$ decay	Decays as $e^{-\alpha r^2}$ , too fast; needs diffuse-function augmentation
One-center (atomic) integrals	Simple, closed-form	Simple, closed-form
Two/three/four-center molecular integrals	Require elliptic coordinates or slow quadrature	Closed-form via Gaussian Product Theorem + Boys function
Number of functions needed per orbital	Few (1–3 STOs can already do well)	Many primitives, but pre-contracted so only 1 "shell" is chosen by the user
Best use case today	Small, precise atomic/diatomic benchmark work; STO-based DFT codes (e.g. ADF, using Slater basis sets efficiently for periodic/molecular DFT)	Essentially all molecular Hartree-Fock, DFT, post-HF, and correlated wavefunction codes
As $\zeta$ or $\alpha \to \infty$	Function collapses toward a delta function at the nucleus (tight core orbital)	Same qualitative collapse, but even more sharply peaked in a Gaussian sense
As $\zeta$ or $\alpha \to 0$	Function spreads over all space; overlap integrals with other centers grow	Same, but with the characteristic Gaussian (not exponential) fall-off

The chemist behind it

The Slater-type orbital is named for John C. Slater, the American physicist whose 1930 rules for estimating effective nuclear charges and orbital exponents gave chemists a simple recipe for constructing physically sensible atomic orbitals without solving the full self-consistent-field problem by hand. Slater worked at a moment when quantum mechanics had just been formalized and the challenge was translating exact hydrogenic mathematics into workable approximations for multi-electron atoms; his exponents ( $\zeta$ values derived from his shielding rules) remain a staple of introductory quantum chemistry courses to this day. Slater's broader career spanned solid-state band theory, the theory of ferromagnetism, and some of the earliest American work connecting quantum theory to digital electronic-structure computation, making him a natural bridge figure between "quantum chemistry on paper" and "quantum chemistry on a machine." The Gaussian-type orbital's decisive role in computation, by contrast, is credited to S. Francis Boys, a British theoretical chemist who in 1950 proposed replacing Slater-type functions with Gaussian functions specifically to make multicenter molecular integrals analytically tractable. Boys recognized — well before electronic computers were widely available to chemists — that the algebraic closure of Gaussians under multiplication (the Gaussian Product Theorem) would convert intractable multicenter integrals into evaluable closed forms, at the cost of needing more functions to represent each atomic orbital adequately. His original insight was, for roughly a decade, more a mathematical curiosity than a practical tool, since hand calculation still made large Gaussian expansions unwieldy; it was the arrival of digital computers in the 1960s and 1970s that turned Boys's proposal into the near-universal foundation of modern quantum chemistry software.

Origins · 1 minJohn Clarke Slater and how this was discovered

John Clarke Slater

1900-1976 · American physicist, Harvard PhD, later MIT professor; pioneer of quantum theory of solids and molecules

Slater developed the Slater determinant (independently formalizing the antisymmetrized many-electron wavefunction), Slater-type orbitals, and Slater's screening rules for effective nuclear charge, while also building foundational theories of energy bands in solids and the exchange interaction underlying ferromagnetism, bridging atomic quantum mechanics and the emerging fields of solid-state and computational physics.

History

- 1900 — John C. Slater born in Oakland, California.

- 1925 — S. Francis Boys born in Pudsey, England (later pioneering Gaussian-basis quantum chemistry).

- 1930 — Slater publishes his screening-constant rules and introduces Slater-type orbitals for atomic and molecular calculations.

- 1950 — Boys proposes using Gaussian functions instead of Slater functions specifically to make molecular integrals analytically tractable.

- 1960s-1970s — Digital computers make large contracted Gaussian basis sets practical, cementing GTOs as the standard for ab initio molecular electronic structure software.

Why can't we just use more Slater functions and avoid Gaussians entirely?

You can, and some codes do — the Amsterdam Density Functional program builds molecular and periodic DFT calculations directly on Slater-type basis functions. The obstacle is not that STOs are inherently unusable; it is that the analytical multicenter two-electron integrals needed for Hartree-Fock and most wavefunction-based post-HF methods do not have closed forms for STOs on four arbitrary centers. STO-based codes instead rely on numerical integration grids or auxiliary fitting schemes, which are efficient for DFT (where the exchange-correlation term is already numerically integrated) but far less competitive for exact-exchange or correlated wavefunction methods, where Gaussian-based analytic integrals dominate for both speed and accuracy control.

What exactly goes wrong if you try to compute a four-center STO integral analytically?

The core problem is geometric: a single STO's level surfaces are spheres about its own nucleus, so the product of two STOs on different nuclei has level surfaces that are prolate spheroids (footballs) around the line joining the two nuclei, not spheres around any single point. For a two-center integral this is still manageable using prolate spheroidal coordinates and tabulated auxiliary functions. But a general two-electron integral $(\mu\nu|\lambda\sigma)$ can involve up to four distinct nuclei simultaneously, and there is no coordinate system that simultaneously spherizes four independent spheroidal constraints. The integral must fall back on multipole expansions or six-dimensional numerical quadrature, both of which scale poorly as molecules grow.

Does using Gaussian functions introduce an error into every quantum chemistry calculation?

Yes, but it is a controllable, quantifiable, and largely removable error, distinct from the theory's own approximations (like the mean-field nature of Hartree-Fock). Because a linear combination of Gaussians can approximate an STO-like shape arbitrarily well as more primitives are added (STO-3G, STO-6G, and beyond), the basis-set error shrinks systematically and predictably with basis size. In practice, chemists quote results at a specific basis-set level (e.g., cc-pVTZ) precisely because the remaining Gaussian-representation error at that level is known to be small and roughly consistent across similar molecules, allowing meaningful comparisons between calculations.

What is the Boys function, physically, and why does it appear everywhere in Gaussian integral codes?

The Boys function $F_0(T) = \int_0^1 e^{-Tt^2}dt$ arises when the Coulomb operator $1/r_{12}$ is rewritten using its Fourier/Gaussian integral representation and then combined with the Gaussian Product Theorem's collapse of multicenter charge distributions into single Gaussian charge clouds. Physically, $F_0(T)$ encodes the Coulomb repulsion energy between two spherical Gaussian charge distributions as a function of a single dimensionless combination of their combined exponent and separation. Because every electron repulsion integral between $s$ -type Gaussian shells reduces to evaluating this one function (and its higher-order analogues $F_n(T)$ for higher angular momentum), efficient, numerically stable algorithms for $F_n(T)$ are a foundational piece of every Gaussian-basis quantum chemistry program.

Is there any situation where Slater-type orbitals genuinely give a better final chemical answer than Gaussians, at fixed cost?

At genuinely fixed number of basis functions per atom, an STO basis typically outperforms a GTO basis of the same size, because a single STO already has the correct cusp and tail that a single GTO can only approximate. This is precisely why STO-3G contracts three Gaussians just to mimic one useful STO. In practice, "at fixed cost" is the wrong comparison, because GTO integrals are so much faster to compute; a Gaussian basis with several times as many primitives per atom, but fully analytic integrals, ultimately reaches a target accuracy faster in wall-clock time than an STO basis relying on slow numerical multicenter integration for a polyatomic molecule.

How does the choice between STOs and GTOs connect to the Roothaan-Hall equations?

The Roothaan-Hall equations turn the Hartree-Fock integro-differential equation into a matrix eigenvalue problem, $\mathbf{FC}=\mathbf{SC}\varepsilon$ , but every entry of the Fock matrix $\mathbf F$ and overlap matrix $\mathbf S$ is itself an integral over the chosen basis functions. If those basis functions are STOs, the two-electron part of $\mathbf F$ requires the difficult multicenter integrals discussed above; if they are Gaussians, those same entries reduce to fast closed-form expressions built from the Gaussian Product Theorem and the Boys function. The practicality of actually solving the Roothaan-Hall equations for a real molecule, rather than merely writing them down, is therefore inseparable from the choice of basis function type.