Configuration Interaction

A single Slater determinant, however carefully optimized, cannot describe what two electrons do when they nearly collide. It gets the average repulsion right but not the instantaneous dodge-and-weave. This is the electron correlation that lowers the true energy below the Hartree–Fock value. Configuration interaction (CI) is the most conceptually direct fix. It builds a bigger, smarter trial wavefunction by mixing the Hartree–Fock determinant with determinants describing electrons promoted into virtual orbitals. A variational calculation then decides how much of each to include. The idea sounds almost too simple to be a real quantum chemistry method — just diagonalize a bigger matrix. That simplicity is exactly why CI is the natural first stop after Hartree–Fock. It converts the many-body correlation problem into ordinary linear algebra. Build the Hamiltonian matrix in a basis of determinants, find its lowest eigenvalue, and read off the eigenvector as the mixing coefficients. Full CI, done in a complete basis, is exact nonrelativistic quantum chemistry. It is the reference against which every cheaper correlation method is judged. But full CI is also a combinatorial catastrophe. The number of determinants grows factorially with the number of electrons and orbitals. So real calculations must truncate the excitation level — singles only (CIS), singles and doubles (CISD), and so on. This truncation buys tractability at a brutal cost: truncated CI is not size-consistent. The energy of two infinitely separated helium atoms computed by CISD is not twice the energy of one helium atom computed by CISD. That failure, first exposed starkly in the 1970s, reshaped the entire field. It motivated the search for methods — coupled cluster and perturbation theory — that get correlation energy without breaking additivity. This article builds CI from the ground up. It covers how the excited determinants arise, and why Brillouin's theorem kills CIS as an energy method. It explains how the secular equations are solved, and exactly how and why size-consistency fails. A fully worked numerical demonstration follows, using a two-electron toy system where the arithmetic can be checked by hand.

Key takeaways
  • CI expands the wavefunction as a linear combination of the HF determinant plus singly, doubly, …, N-tuply excited determinants, with coefficients found variationally.
  • Brillouin's theorem shows singly-excited determinants do not couple to the HF reference through H^\hat{H}, so CIS gives zero correlation energy correction to the ground state (but describes excited states).
  • CISD (singles + doubles) recovers most dynamic correlation for small systems but misses higher excitations that matter for larger systems.
  • Full CI in a given one-electron basis is the exact solution of the electronic Schrödinger equation within that basis — the ultimate benchmark.
  • Truncated CI (CID, CISD, …) is not size-consistent: ECISD(A⋯B)≠ECISD(A)+ECISD(B)E_{\text{CISD}}(A\cdots B) \neq E_{\text{CISD}}(A) + E_{\text{CISD}}(B) for noninteracting fragments, because disconnected quadruple excitations are missing.
  • Key idea: CI turns "find the correlated wavefunction" into "diagonalize a Hamiltonian matrix in a basis of Slater determinants," and the truncation needed to make that matrix small is precisely what breaks additivity of the energy.

Where this lives today

Full CI, though exponentially expensive, is still the gold-standard benchmark used to calibrate every approximate correlation method — quantum chemists routinely run full CI on tiny model systems (a handful of electrons in a modest basis) purely to certify that a new cheaper method (a new density functional, a new coupled-cluster truncation, a new machine-learned potential) gets the right answer where the exact one is computable. CISD and its multireference cousin MRCI (multireference CI, often with a Davidson size-consistency correction added ad hoc) remain the workhorses for excited-state spectroscopy and photochemistry — computing UV-Vis absorption spectra of organic dyes, conical intersections in photoisomerization (the mechanism behind vision, via retinal photochemistry), and transition-metal d–d and charge-transfer transitions where multiple electronic configurations are intrinsically close in energy. Selected/adaptive CI algorithms (CIPSI, Heat-bath CI, semistochastic FCI-QMC) are a modern renaissance of the CI idea: rather than truncating by excitation level, they truncate by importance, keeping only the determinants that actually matter energetically, and have pushed near-exact CI-quality benchmarks to systems with dozens of correlated electrons — directly informing where cheaper production methods like density functional theory and coupled cluster theory (see Moller-Plesset Perturbation Theory (MP2) and Coupled Cluster Theory for the size-consistent alternatives this article's failure motivates) can be trusted or must be corrected.

Ask the tutor: Why does representing the wavefunction as an exponential of excitation operators (as in coupled cluster theory) automatically restore size-consistency where a truncated linear CI expansion cannot?

The physical picture

Start from the Hartree–Fock ground state ∣Ψ0⟩|\Psi_0\rangle, a single Slater determinant built from the lowest NN spin orbitals out of a full one-electron basis of KK spin orbitals (so there are K−NK-N empty, "virtual" orbitals). Hartree–Fock is the best single determinant, found by the variational optimization discussed in Electron Correlation and the Correlation Energy and built on the antisymmetry requirement of Antisymmetry and the Slater Determinant. But the true ground state is not a single determinant. Configuration interaction restores flexibility by writing

∣ΨCI⟩=c0∣Ψ0⟩+∑arcar∣Ψar⟩+∑a<b, r<scabrs∣Ψabrs⟩+⋯|\Psi_{\text{CI}}\rangle = c_0|\Psi_0\rangle + \sum_{ar} c_a^r |\Psi_a^r\rangle + \sum_{a<b,\,r<s} c_{ab}^{rs} |\Psi_{ab}^{rs}\rangle + \cdots

where ∣Ψar⟩|\Psi_a^r\rangle denotes the determinant obtained from ∣Ψ0⟩|\Psi_0\rangle by replacing occupied spin orbital aa with virtual spin orbital rr (a "single excitation"), ∣Ψabrs⟩|\Psi_{ab}^{rs}\rangle replaces two occupied orbitals with two virtuals ("double excitation"), and so on up to NN-fold excitations. Each determinant is itself a full antisymmetrized product of spin orbitals — nothing but book-keeping on top of Slater's construction — but the combination is no longer a single product: it is a genuine many-body superposition. Physically, the double excitations are the workhorses. They describe the instantaneous correlation of electron pairs — when electron 1 is close to electron 2 in the HF description, exciting both into virtuals with the right phase lets the wavefunction build in a "hole" around each electron, the Coulomb correlation hole. Singles, by contrast, mostly relax orbitals (they mix in a little orbital rotation) and contribute almost nothing to the ground-state correlation energy directly, for a reason made precise by Brillouin's theorem below. Triples and quadruples matter increasingly as the system gets bigger, because they can be built approximately as products of independent double excitations on separate parts of the molecule — and it is exactly this product structure that truncated CI cannot represent, since CISD stops at pairs and cannot see the simultaneous double-double excitation on two well-separated fragments.

Figure
Determinant hierarchy: HF reference and excitations into virtual orbitals
Excitation levels above the Hartree–Fock determinant Ψ₀ (HF) occ. a,b virt. r,s Ψ ᵣₐ (single) Ψ ᵣˢₐᵦ (double) Ψ (quad.) CISD keeps only |Ψ₀⟩, singles, doubles. A disconnected "double-on-A × double-on-B" quadruple (far right) is invisible to CISD — this is the microscopic origin of the size-consistency failure.

The variational principle guarantees that adding more determinants can only lower (or leave unchanged) the ground-state energy, since CI is nothing but a linear variational calculation in an expanded basis of many-electron functions. The lowest eigenvalue of the CI Hamiltonian matrix is therefore always an upper bound to the true ground-state energy, and it approaches the exact answer (within the one-electron basis) as more excitation classes are included — reaching it exactly at full CI.

Figure
Energy lowering as excitation levels are added, schematically
Variational energy vs. excitation level included excitation level included → E (hartree, schematic) E(exact, full CI) HF +S +D +T +Q Full CI

The core ideas

The CI expansion and secular equations. Writing the trial wavefunction as a linear combination of all NN-electron Slater determinants {∣ΨI⟩}\{|\Psi_I\rangle\} built from a fixed one-electron basis,

∣ΨCI⟩=∑IcI∣ΨI⟩,|\Psi_{\text{CI}}\rangle = \sum_I c_I |\Psi_I\rangle ,

where: I=0I=0 labels the HF reference and I>0I>0 labels singly, doubly, … excited determinants, and cIc_I are variational coefficients — minimizing the Rayleigh quotient ⟨ΨCI∣H^∣ΨCI⟩/⟨ΨCI∣ΨCI⟩\langle \Psi_{\text{CI}}|\hat H|\Psi_{\text{CI}}\rangle / \langle \Psi_{\text{CI}}|\Psi_{\text{CI}}\rangle subject to normalization gives the matrix eigenvalue problem

Hc=E c,HIJ=⟨ΨI∣H^∣ΨJ⟩.\mathbf{H}\mathbf{c} = E\,\mathbf{c}, \qquad H_{IJ} = \langle \Psi_I | \hat H | \Psi_J \rangle .

where: H\mathbf{H} is the CI (secular) matrix, EE its eigenvalues (ground and excited electronic states simultaneously), c\mathbf c the eigenvector of expansion coefficients. This is the exact many-body analogue of the linear variation method used for simple trial functions — only now the "basis functions" are entire Slater determinants. Reading: diagonalizing H\mathbf H once yields the entire spectrum: the lowest eigenvalue is the correlated ground state, and higher eigenvalues are approximate excited states — CI is simultaneously a ground-state and excited-state method, unlike Hartree–Fock.

Slater–Condon rules. The matrix elements HIJH_{IJ} are never computed by brute-force integration over NN-electron coordinates; they follow from compact rules depending on how many spin orbitals differ between ΨI\Psi_I and ΨJ\Psi_J:

⟨Ψ0∣H^∣Ψar⟩=⟨a∣f^∣r⟩=far=0(Brillouin’s theorem, HF orbitals),\langle \Psi_0|\hat H|\Psi_a^r\rangle = \langle a|\hat f|r\rangle = f_{ar} = 0 \quad (\text{Brillouin's theorem, HF orbitals}),
⟨Ψ0∣H^∣Ψabrs⟩=⟨ab∣∣rs⟩,\langle \Psi_0|\hat H|\Psi_{ab}^{rs}\rangle = \langle ab||rs\rangle ,
⟨Ψ0∣H^∣Ψabcrst⟩=0.\langle \Psi_0|\hat H|\Psi_{abc}^{rst}\rangle = 0 .

where: f^\hat f is the Fock operator whose eigenfunctions are the HF orbitals, farf_{ar} its off-diagonal occupied–virtual element (zero at self-consistency), and ⟨ab∣∣rs⟩=⟨ab∣rs⟩−⟨ab∣sr⟩\langle ab||rs\rangle = \langle ab|rs\rangle - \langle ab|sr\rangle is an antisymmetrized two-electron integral in physicists' notation. Reading: Brillouin's theorem is the reason CIS (which only adds single excitations to Ψ0\Psi_0) gives a ground-state energy that is identical to Hartree–Fock through first order — singles do not couple to Ψ0\Psi_0 directly. The Hamiltonian is a two-body operator, so Ψ0\Psi_0 couples directly only to doubles; triples and higher couple to Ψ0\Psi_0 only indirectly, through their coupling to doubles and singles inside the CI matrix. This is what makes doubles the dominant correlating configurations and is the microscopic reason CISD is the traditional workhorse truncation.

Size-consistency (formal statement). For two well-separated, noninteracting subsystems AA and BB, exactness requires

E(A⋯B)=E(A)+E(B).E(A\cdots B) = E(A) + E(B).

where: A⋯BA\cdots B denotes the supersystem at a separation large enough that all inter-fragment interactions vanish. Full CI satisfies this identically; any truncated CI (CID, CISD, CISDT, …) generally does not, because a simultaneous double excitation on AA and a double excitation on BB is a quadruple excitation relative to the supersystem reference, and CISD excludes quadruples.

Derivations

Derivation ABrillouin's theorem: why singles do not mix into the HF ground state

Goal. Show that ⟨Ψ0∣H^∣Ψar⟩=0\langle \Psi_0 | \hat H | \Psi_a^r \rangle = 0 when Ψ0\Psi_0 is built from canonical Hartree–Fock orbitals, and understand the consequence for CIS.

Step 1 — Write the electronic Hamiltonian in second-quantized/matrix-element form. For NN electrons in spin orbitals,

H^=∑ih^(i)+∑i<j1rij,\hat H = \sum_i \hat h(i) + \sum_{i<j} \frac{1}{r_{ij}} ,

where h^(i)\hat h(i) is the one-electron core Hamiltonian (kinetic energy plus nuclear attraction) for electron ii. Using the Slater–Condon rule for determinants differing by exactly one spin orbital (a→ra \to r),

⟨Ψ0∣H^∣Ψar⟩=⟨a∣h^∣r⟩+∑b (occ)⟨ab∣∣rb⟩.\langle \Psi_0|\hat H|\Psi_a^r\rangle = \langle a|\hat h|r\rangle + \sum_{b\,(\text{occ})} \langle ab||rb\rangle .

Step 2 — Recognize the Fock matrix element. The Fock operator acting on spin orbital aa is defined, for an NN-electron closed-shell or spin-orbital HF reference, as

f^(1)=h^(1)+∑b (occ)[J^b(1)−K^b(1)],\hat f(1) = \hat h(1) + \sum_{b\,(\text{occ})} \big[\hat J_b(1) - \hat K_b(1)\big] ,

so that

⟨a∣f^∣r⟩=⟨a∣h^∣r⟩+∑b⟨ab∣∣rb⟩.\langle a|\hat f|r\rangle = \langle a|\hat h|r\rangle + \sum_{b} \langle ab||rb\rangle .

Comparing with Step 1, the right-hand side is exactly ⟨a∣f^∣r⟩=far\langle a | \hat f | r \rangle = f_{ar}:

⟨Ψ0∣H^∣Ψar⟩=far.\langle \Psi_0|\hat H|\Psi_a^r\rangle = f_{ar}.

Step 3 — Use the defining property of canonical HF orbitals. Canonical Hartree–Fock orbitals are chosen (by a unitary rotation among occupied orbitals and separately among virtuals) to diagonalize the Fock matrix within each block, and the self-consistent-field condition additionally forces the occupied–virtual block to vanish:

far=0for all occupied a, virtual r.f_{ar} = 0 \quad \text{for all occupied } a,\ \text{virtual } r .

This is a direct consequence of the HF variational condition δEHF/δ{orbital rotations}=0\delta E_{\text{HF}}/\delta \{\text{orbital rotations}\} = 0 applied specifically to occupied–virtual rotations.

Step 4 — Conclude. Therefore

⟨Ψ0∣H^∣Ψar⟩=0\langle \Psi_0 | \hat H | \Psi_a^r \rangle = 0

for every single excitation, when the orbitals used are the converged canonical HF set. Consequently, in the CI secular matrix, the entire first row/column connecting Ψ0\Psi_0 to all singles is exactly zero. By first-order perturbation theory this means singles contribute zero first-order correction to the HF ground-state energy and, by Löwdin partitioning, their leading effect on E0E_0 enters only at second order through mixing mediated by doubles. This is why a CI truncated at singles only (CIS) reproduces EHFE_{\text{HF}} for the ground state and is used almost exclusively for excited states, where the singles do mix with each other and with doubles.

Go deeper · Derivation B · 2 minSize-inconsistency of truncated CI: the two-noninteracting-helium argument
Derivation BSize-inconsistency of truncated CI: the two-noninteracting-helium argument

Goal. Prove explicitly that CID (doubles-only CI, the cleanest case) applied to a supersystem of two infinitely separated, identical two-electron atoms does not give twice the monomer CID energy — and identify exactly which configurations are missing.

Step 1 — Set up the monomer problem. Let atom AA be a two-electron system (e.g., a minimal-basis He atom) with HF ground state ∣Ψ0A⟩|\Psi_0^A\rangle built from occupied orbital aa and some virtual orbitals. In a minimal basis with one occupied and one virtual spin-orbital pair (a 2-orbital model, "full CI on He" toy), the monomer CID wavefunction is

∣ΨA⟩=∣Ψ0A⟩+cA ∣Ψaa′rr′⟩,|\Psi^A\rangle = |\Psi_0^A\rangle + c^A\,|\Psi_{aa'}^{rr'}\rangle ,

where Ψaa′rr′\Psi_{aa'}^{rr'} is the double excitation of both electrons of atom AA from the occupied pair into the virtual pair. Solving the 2×22\times2 CID secular equation for the monomer gives ground-state energy ECIDA=E0A+ΔAE^A_{\text{CID}} = E_0^A + \Delta^A for some correlation lowering ΔA<0\Delta^A<0, exact in this minimal model because with only one occupied spatial orbital, CID is full CI for the monomer.

Step 2 — Build the supersystem reference at infinite separation. Put atom BB, an identical copy, infinitely far from AA so there is no interaction: H^AB=H^A+H^B\hat H_{AB} = \hat H_A + \hat H_B. The exact ground state factorizes,

Ψexact(A,B)=ΨexactA×ΨexactB,\Psi_{\text{exact}}(A,B) = \Psi^A_{\text{exact}} \times \Psi^B_{\text{exact}},

and the exact total energy is additive: Eexact=EexactA+EexactBE_{\text{exact}} = E^A_{\text{exact}} + E^B_{\text{exact}}. Expand each factor as reference plus double excitation, ΨexactA≈Ψ0A+cAΨaa′rr′(A)\Psi^A_{\text{exact}} \approx \Psi_0^A + c^A \Psi_{aa'}^{rr'}(A) and similarly for BB. Multiplying out,

Ψexact(A,B)≈Ψ0AΨ0B+cA Ψaa′rr′(A)Ψ0B+cB Ψ0A Ψaa′rr′(B)+cAcB Ψaa′rr′(A) Ψaa′rr′(B).\Psi_{\text{exact}}(A,B) \approx \Psi_0^A\Psi_0^B + c^A\,\Psi_{aa'}^{rr'}(A)\Psi_0^B + c^B\,\Psi_0^A\,\Psi_{aa'}^{rr'}(B) + c^A c^B\, \Psi_{aa'}^{rr'}(A)\,\Psi_{aa'}^{rr'}(B).

Step 3 — Classify each term by excitation level relative to the supersystem HF reference Ψ0AΨ0B\Psi_0^A\Psi_0^B. The first term is the reference itself. The second and third terms are double excitations (two electrons on AA promoted, BB untouched, or vice versa) — CID keeps these. The fourth term, cAcB Ψaa′rr′(A)Ψaa′rr′(B)c^Ac^B\,\Psi_{aa'}^{rr'}(A)\Psi_{aa'}^{rr'}(B), has four spin orbitals changed relative to the supersystem reference (two on AA, two on BB simultaneously): it is a quadruple excitation. Crucially, it is a disconnected quadruple — simply the product of two independent doubles on non-interacting fragments — and its coefficient factorizes exactly as cAcBc^Ac^B, the product of the monomer amplitudes. There is nothing exotic about it; it is required by simple algebra (multiplying out a product of two sums).

Step 4 — Show CID (or CISD) supersystem excludes exactly this term. CID for the supersystem is defined as Ψ0AΨ0B\Psi_0^A\Psi_0^B plus all double excitations relative to that reference — but by definition it stops at doubles, so the quadruple term cAcB Ψaa′rr′(A)Ψaa′rr′(B)c^Ac^B\,\Psi^{rr'}_{aa'}(A)\Psi^{rr'}_{aa'}(B) is not included in the CID variational space at all. Solving the supersystem CID secular equations with only reference + doubles necessarily gives an energy

ECID(A⋯B)≠ECID(A)+ECID(B),E_{\text{CID}}(A\cdots B) \neq E_{\text{CID}}(A) + E_{\text{CID}}(B),

and in fact (see the worked example below) gives an energy higher than the additive sum — the variational calculation is starved of a term the exact wavefunction needs, and by the variational principle any energy computed with a wavefunction missing a needed component cannot be lower than one computed with it present. Since the missing disconnected quadruple is exactly the product term forced by simple factorization, the size-consistency error grows with the number of noninteracting units — it is not a small correction, it is a qualitative failure that gets worse as the "molecule" (or cluster of identical fragments) gets larger. This single derivation is why truncated CI was abandoned as a general-purpose method for large systems in favor of size-consistent alternatives such as coupled-cluster theory.

Worked examples

Example 1 — Foundational: building and diagonalizing a minimal CI matrix by hand (two-electron, two-orbital model). Take a toy two-electron "atom" with one occupied spatial orbital ϕ1\phi_1 (HF occupied, energy ε1\varepsilon_1) and one virtual spatial orbital ϕ2\phi_2 (energy ε2\varepsilon_2), doubly occupied/empty in the closed-shell sense (spatial CI, singlet-adapted, so we need only the reference ∣Ψ0⟩=∣ϕ1ϕˉ1∣|\Psi_0\rangle = |\phi_1\bar\phi_1| and the doubly-excited singlet ∣Ψ2⟩=∣ϕ2ϕˉ2∣|\Psi_2\rangle = |\phi_2\bar\phi_2|; the singly-excited configuration vanishes by symmetry/Brillouin's theorem for a closed shell). Use these representative numbers (not from any particular real atom, chosen asymmetric and dimensionless in hartree for clarity):

ε1=−1.4,ε2=−0.30,K12=⟨12∣12⟩=0.62,⟨11∣22⟩=0.35.\varepsilon_1 = -1.4,\qquad \varepsilon_2 = -0.30,\qquad K_{12} = \langle 12|12\rangle = 0.62,\qquad \langle 11|22\rangle = 0.35 .

The 2×22\times2 CI matrix in the basis {Ψ0,Ψ2}\{\Psi_0, \Psi_2\} (singlet-coupled double excitation) is standard textbook algebra:

H00=2ε1+J11,H22=2ε2+J22,H02=⟨11∣22⟩.H_{00} = 2\varepsilon_1 + J_{11}, \qquad H_{22} = 2\varepsilon_2 + J_{22}, \qquad H_{02} = \langle 11|22\rangle .

Take J11=0.68J_{11} = 0.68 (Coulomb self-repulsion of ϕ1\phi_1) and J22=0.42J_{22} = 0.42. Then

H00=2(−1.4)+0.68=−2.12,H22=2(−0.30)+0.42=−0.18.H_{00} = 2(-1.4) + 0.68 = -2.12, \qquad H_{22} = 2(-0.30) + 0.42 = -0.18 .

The secular determinant is

∣−2.12−E0.350.35−0.18−E∣=0  ⟹  (E+2.12)(E+0.18)−0.1225=0.\begin{vmatrix} -2.12 - E & 0.35 \\ 0.35 & -0.18-E \end{vmatrix} = 0 \;\Longrightarrow\; (E+2.12)(E+0.18) - 0.1225 = 0.

Expanding:

E2+2.30E+0.3816−0.1225=0  ⟹  E2+2.30E+0.2591=0.E^2 + 2.30E + 0.3816 - 0.1225 = 0 \;\Longrightarrow\; E^2 + 2.30E + 0.2591 = 0.

The quadratic formula gives

E=−2.30±2.302−4(0.2591)2=−2.30±5.29−1.03642=−2.30±4.25362.E = \frac{-2.30 \pm \sqrt{2.30^2 - 4(0.2591)}}{2} = \frac{-2.30 \pm \sqrt{5.29 - 1.0364}}{2} = \frac{-2.30 \pm \sqrt{4.2536}}{2}.
4.2536=2.0624,E−=−2.30−2.06242=−2.1812,E+=−2.30+2.06242=−0.1188.\sqrt{4.2536} = 2.0624, \qquad E_- = \frac{-2.30 - 2.0624}{2} = -2.1812, \qquad E_+ = \frac{-2.30+2.0624}{2} = -0.1188.

The lower root, ECI=−2.1812E_{\text{CI}} = -2.1812 hartree, is the (full) CI ground-state energy — lower than H00=−2.12H_{00}=-2.12, the Hartree–Fock energy, by

Ecorr=ECI−H00=−2.1812−(−2.12)=−0.0612 hartree,E_{\text{corr}} = E_{\text{CI}} - H_{00} = -2.1812 - (-2.12) = -0.0612\ \text{hartree},

the correlation energy recovered by mixing in the double excitation. Because this model has only one occupied and one virtual spatial orbital, this doubles-only CI is full CI for this system — an explicit, hand-checkable instance of the variational lowering promised by the theory. As a limiting-case check: if the coupling ⟨11∣22⟩→0\langle 11|22\rangle \to 0 (no correlation term), the secular equation decouples and ECI→H00=−2.12E_{\text{CI}} \to H_{00} = -2.12 exactly, recovering the Hartree–Fock result with zero correlation energy, as it must.

Example 2 — Applied: CISD vs. CID vs. full CI energy ordering for a stretched bond (qualitative but numerically grounded). Consider H2_2 described in a minimal basis (this is the smallest system where singles, doubles, and the size-consistency question can all be illustrated together, and where full CI is exactly CISD — there are no triples or quadruples possible with only two electrons). In a minimal (STS-3G–like, two-orbital) basis, H2_2 has one bonding MO σg\sigma_g (occupied) and one antibonding MO σu∗\sigma_u^* (virtual); the only excitations available are the double excitation σg2→σu∗2\sigma_g^2 \to \sigma_u^{*2} (the single excitation vanishes by symmetry, since σg→σu∗\sigma_g \to \sigma_u^* changes overall spatial symmetry of a closed-shell singlet incorrectly — it is a triplet, not coupling to the singlet ground state). At the equilibrium bond length (R≈0.74R \approx 0.74 Å), representative minimal-basis numbers are

EHF≈−1.1167 Eh,EFCI≈−1.1367 Eh,Ecorr≈−0.020 Eh.E_{\text{HF}} \approx -1.1167\ E_h, \qquad E_{\text{FCI}} \approx -1.1367\ E_h, \qquad E_{\text{corr}} \approx -0.020\ E_h.

At a stretched bond (R≈3.0R \approx 3.0 Å, well past the point where the ionic HF description qualitatively fails), the double excitation coefficient c2c_2 becomes large (approaching c2→−1/2c_2 \to -1/\sqrt2 relative to c0→1/2c_0 \to 1/\sqrt2 as R→∞R\to\infty, since the true wavefunction must become an equal mix of covalent and ionic-canceling configurations to dissociate correctly into two neutral H atoms). Representative stretched-bond numbers:

EHF(3.0 A˚)≈−0.70 Eh,EFCI(3.0 A˚)≈−1.00 Eh,E_{\text{HF}}(3.0\,\text{Å}) \approx -0.70\ E_h, \qquad E_{\text{FCI}}(3.0\,\text{Å}) \approx -1.00\ E_h,

so the correlation energy grows to roughly −0.30 Eh-0.30\ E_h — nearly fifteen times its equilibrium value — because Hartree–Fock's single-determinant description becomes qualitatively wrong at dissociation (it retains 50% ionic H+^+H−^- character that has no business surviving to infinite separation), while full CI corrects this exactly since the double excitation is precisely the configuration needed to cancel the spurious ionic term. This is why static correlation (large coefficient on one specific excited determinant, here the double) becomes dominant at stretched geometries, in contrast to the dynamic correlation (many small coefficients) that dominates near equilibrium — CI captures both, but a truncated CI or a perturbative method tuned for the equilibrium (dynamic-correlation) regime often fails badly in the stretched (static-correlation) regime. Because H2_2 has exactly two electrons, CISD = full CI here and there is no size-consistency error for the isolated molecule — the size-consistency problem appears only once we ask about two separated H2_2 molecules simultaneously, which is Example 3.

Example 3 — Competition-level: quantifying the size-consistency error for two identical noninteracting systems. Using the results of Example 1, suppose we now have two copies of that toy two-electron atom, infinitely far apart, each individually described exactly by the 2×22\times2 CI of Example 1 (ECIDmonomer=−2.1812 EhE^{\text{monomer}}_{\text{CID}} = -2.1812\ E_h each). The exact (and, since each monomer's "CID" is already full CI, also the exact full-CI) energy of the noninteracting pair is additive:

Eexact(A+B)=ECIDA+ECIDB=2(−2.1812)=−4.3624 Eh.E_{\text{exact}}(A+B) = E^A_{\text{CID}} + E^B_{\text{CID}} = 2(-2.1812) = -4.3624\ E_h.

Now perform CID (doubles only, relative to the supersystem reference Ψ0AΨ0B\Psi_0^A\Psi_0^B) directly on the 2-fragment supersystem, allowing only single-fragment double excitations Ψaa′rr′(A)Ψ0B\Psi_{aa'}^{rr'}(A)\Psi_0^B and Ψ0AΨaa′rr′(B)\Psi_0^A\Psi_{aa'}^{rr'}(B) — not the simultaneous double-double (quadruple) term. By symmetry (identical fragments) the supersystem CID secular problem is a 3×33\times3 matrix in the basis {Ψ0AΨ0B, Ψ2AΨ0B, Ψ0AΨ2B}\{\Psi_0^A\Psi_0^B,\ \Psi_2^A\Psi_0^B,\ \Psi_0^A\Psi_2^B\}:

Hsuper=(2H00H02H02H02H00+H220H020H00+H22)=(−4.240.350.350.35−2.3000.350−2.30).\mathbf H_{\text{super}} = \begin{pmatrix} 2H_{00} & H_{02} & H_{02} \\ H_{02} & H_{00}+H_{22} & 0 \\ H_{02} & 0 & H_{00}+H_{22} \end{pmatrix} = \begin{pmatrix} -4.24 & 0.35 & 0.35\\ 0.35 & -2.30 & 0\\ 0.35 & 0 & -2.30 \end{pmatrix}.

(The off-diagonal zero between the two single-fragment doubles holds because H^A\hat H_A and H^B\hat H_B act on disjoint electrons — this block structure is itself only exactly this simple in the noninteracting limit.) By symmetry, one eigenvector is the antisymmetric combination (Ψ2AΨ0B−Ψ0AΨ2B)/2(\Psi_2^A\Psi_0^B - \Psi_0^A\Psi_2^B)/\sqrt2 with eigenvalue exactly H00+H22=−2.30 EhH_{00}+H_{22} = -2.30\ E_h (it does not mix with the reference, since the two off-diagonal 0.350.35's cancel). The relevant symmetric 2×22\times2 block, in the basis {Ψ0AΨ0B, (Ψ2AΨ0B+Ψ0AΨ2B)/2}\{\Psi_0^A\Psi_0^B,\ (\Psi_2^A\Psi_0^B+\Psi_0^A\Psi_2^B)/\sqrt2\}, has off-diagonal element 2×0.35=0.4950\sqrt2 \times 0.35 = 0.4950:

(−4.240.49500.4950−2.30).\begin{pmatrix} -4.24 & 0.4950 \\ 0.4950 & -2.30 \end{pmatrix}.

Solving:

(E+4.24)(E+2.30)−0.2450=0  ⟹  E2+6.54E+9.752−0.2450=0  ⟹  E2+6.54E+9.507=0.(E+4.24)(E+2.30) - 0.2450 = 0 \;\Longrightarrow\; E^2 + 6.54E + 9.752 - 0.2450 = 0 \;\Longrightarrow\; E^2+6.54E+9.507=0.
E=−6.54±6.542−4(9.507)2=−6.54±42.7716−38.0282=−6.54±4.74362.E = \frac{-6.54 \pm \sqrt{6.54^2 - 4(9.507)}}{2} = \frac{-6.54\pm\sqrt{42.7716-38.028}}{2} = \frac{-6.54\pm\sqrt{4.7436}}{2}.
4.7436=2.1780,E−=−6.54−2.17802=−4.359.\sqrt{4.7436}=2.1780, \qquad E_- = \frac{-6.54-2.1780}{2} = -4.359.

So ECIDsupersystem=−4.359 EhE_{\text{CID}}^{\text{supersystem}} = -4.359\ E_h, compared with the exact additive value −4.3624 Eh-4.3624\ E_h found above. The size-consistency error is

ΔSC=ECIDsupersystem−2ECIDmonomer=−4.359−(−4.3624)=+0.0034 Eh≈+2.1 kcal/mol,\Delta_{\text{SC}} = E_{\text{CID}}^{\text{supersystem}} - 2E^{\text{monomer}}_{\text{CID}} = -4.359 - (-4.3624) = +0.0034\ E_h \approx +2.1\ \text{kcal/mol},

a genuine, nonzero, positive error (CID for the pair recovers less than twice the per-monomer correlation energy) purely because the disconnected quadruple excitation Ψ2AΨ2B\Psi_2^A\Psi_2^B was excluded from the variational space — exactly the mechanism identified in the derivation above. As a limiting-case check, if the two monomers were allowed to interact (HAB≠0H_{AB}\neq0, the two-monomer basis would also need genuine four-electron double excitations spanning both fragments together, which are included at the CID/CISD level for a "real" bonded 4-electron molecule); it is specifically the noninteracting, well-separated limit that isolates the pure size-consistency artifact from ordinary intermolecular correlation, and the error shown here — small in this toy model but structurally the same mechanism — grows linearly with the number of identical noninteracting units in general, eventually dominating for large clusters or polymers.

Limiting cases, scope & common traps

Regime / limitWhat happensWhy
Full CI, any systemExact nonrelativistic energy within the one-electron basisAll excitation levels present; no missing determinant, so no truncation error at all
CIS (singles only)Ground-state energy = EHFE_{\text{HF}} exactlyBrillouin's theorem: ⟨Ψ0∣H^∣Ψar⟩=0\langle\Psi_0|\hat H|\Psi_a^r\rangle=0 for canonical HF orbitals
CID / CISDRecovers most dynamic correlation for small, compact moleculesDoubles couple directly to Ψ0\Psi_0; but disconnected quadruples on separated fragments are missing
Two or more noninteracting fragmentsECISD(A+B)≠ECISD(A)+ECISD(B)E_{\text{CISD}}(A{+}B) \neq E_{\text{CISD}}(A)+E_{\text{CISD}}(B)Missing disconnected higher excitations (Example 3)
Bond stretching / dissociationStatic correlation dominates; HF and low-order truncated CI both degrade, but CI (with enough excitation levels) still improves qualitativelyNear-degenerate configurations (e.g. the doubly excited determinant) become comparably weighted to Ψ0\Psi_0
N→∞N\to\infty (large molecule, extended system)Truncated CI errors grow catastrophically; size-consistency failure compoundsMore independent "fragments" each contribute missing disconnected excitations
Two-electron systems (e.g. H2_2, He)CISD = Full CI automaticallyOnly doubles (and, if symmetry allows, no singles) exist; no triples/quadruples possible with two electrons

The chemist behind it

Origins · 1 minPer-Olov Löwdin and how this was discovered
PL
Per-Olov Löwdin
1916–2000 · Quantum Theory Project, Uppsala & University of Florida

Löwdin was one of the principal architects of postwar quantum chemistry as a mathematically rigorous discipline, formalizing natural orbitals, symmetric orthogonalization, and the partitioning technique used to analyze CI matrices and truncation errors. He founded the Sanibel Symposia and the Quantum Theory Project, institutions that turned configuration interaction and related many-body methods from isolated calculations into a systematic international research program. His insistence on precise variational bounds and natural-orbital analysis directly shaped how chemists came to understand — and eventually solve — the size-consistency problem of truncated CI.

History
- 1929 — Configuration interaction is anticipated in early atomic spectroscopy: mixing of nearly degenerate atomic configurations (term interaction) explains anomalous spectral line intensities and positions.
- 1955 — Löwdin formalizes natural orbitals and partitioning techniques, giving CI a rigorous variational and perturbative footing.
- 1960s — Full CI and truncated CI calculations become feasible on early digital computers for small molecules (H2_2, LiH, BeH2_2), establishing CISD as the standard correlated method of the era.
- 1974–1977 — Explicit demonstrations (notably by Pople and coworkers, and Bartlett and Purvis) show CISD is not size-consistent, using dissociation and cluster-of-atoms tests; the result reshapes correlation-method development.
- 1980s–1990s — Coupled cluster theory (CCSD, CCSD(T)) and size-consistent multireference CI variants (e.g. Davidson-corrected CISD, MRCI+Q) become the practical response to the size-consistency problem.
Why does CI use Slater determinants rather than simple orbital products as the expansion basis?

Electrons are fermions, so the total wavefunction must be antisymmetric under exchange of any two electrons' coordinates — a requirement built directly into the Slater determinant construction. Using determinants as the building blocks of CI guarantees every term in the expansion, and hence the full CI wavefunction, automatically satisfies the Pauli antisymmetry requirement without any extra bookkeeping. A naive product expansion would need explicit antisymmetrization glued on afterward, and would double-count or omit exchange effects. Determinants are also the natural output of Hartree–Fock, so CI can be built directly on top of an existing, already-optimized HF calculation, reusing its orbitals and orbital energies with no extra machinery.

Why doesn't CIS improve the ground-state energy at all, if it clearly changes the wavefunction?

Brillouin's theorem shows that, using canonical Hartree–Fock orbitals, every singly-excited determinant has exactly zero Hamiltonian matrix element with the reference determinant. In a CI matrix restricted to {Ψ0,singles}\{\Psi_0,\text{singles}\}, the reference row/column coupling to singles is entirely zero, so diagonalizing the block containing Ψ0\Psi_0 returns EHFE_{\text{HF}} unchanged as the lowest eigenvalue coupled to that reference. Singles do mix strongly with each other and can produce meaningful excited-state energies (CIS is a standard cheap method for UV-Vis excitation energies), but for the ground state specifically, singles are a "dead end" without doubles present to mediate any indirect coupling.

Is full CI ever actually used for real molecules, or is it purely a theoretical benchmark?

Full CI's cost scales combinatorially — roughly as the binomial coefficient of orbitals choose electrons for each spin — so exact full CI is limited in practice to small molecules in modest basis sets (systems with perhaps 10–18 correlated electrons in a double-zeta basis, depending on available computer resources). Within that regime it is used as a genuine, not merely theoretical, benchmark: method developers routinely run full CI on molecules like water, methylene, or small hydrides to certify new density functionals, new coupled-cluster truncations, or new machine-learned interatomic potentials against the numerically exact nonrelativistic answer. Selected-CI and stochastic FCI-QMC methods have extended this exact benchmarking capability to considerably larger systems in recent years.

How is the size-consistency problem actually fixed in modern quantum chemistry?

Two main strategies are used. The first is to abandon the truncated-CI linear expansion in favor of an exponential parametrization of the wavefunction — coupled cluster theory — where the cluster operator automatically generates the necessary disconnected higher excitations as products of lower ones, restoring exact size-consistency (and, for the standard truncations like CCSD, full size-extensivity) essentially for free. The second is empirical correction of CISD itself, most commonly the Davidson correction, which estimates the missing quadruple-excitation contribution from the CISD singles/doubles coefficients and adds it back approximately. Coupled cluster has become the dominant modern solution precisely because it fixes the problem structurally rather than patching it after the fact.

Does size-consistency matter for a single, non-dissociating molecule, or only for comparing separated fragments?

It matters even for one molecule whenever you compare energies across different numbers of "independent" correlating units within that molecule — for example, computing a large hydrocarbon's total correlation energy, or comparing a dimer's binding energy against two monomers computed separately (a very common calculation in noncovalent-interaction and supramolecular chemistry). Any calculation that implicitly or explicitly compares energies at different particle numbers, or across different-sized fragments of an extended system, is exposed to the size-consistency error if a truncated CI method is used. For a single fixed molecule computed once, the absolute error is present but is usually absorbed into the general truncation error rather than showing up as a separate, obviously diagnosable artifact.

Why do doubles couple directly to the HF reference but triples do not?

The electronic Hamiltonian contains only one-electron and two-electron operators — it can, at most, simultaneously change the occupation of two spin orbitals in a single matrix element (via the two-electron repulsion term) or one spin orbital (via the one-electron term). A determinant differing from the reference in three spin orbitals (a triple excitation) therefore has an exactly zero direct Hamiltonian matrix element with Ψ0\Psi_0, by the Slater–Condon rules, regardless of the orbital choice — this is a much more general and robust rule than Brillouin's theorem (which needed the extra assumption of canonical HF orbitals). Triples still matter energetically because they mix with doubles and singles inside the full CI matrix, providing an indirect but often significant contribution.

Is CI compatible with density functional theory, or are they alternative, mutually exclusive frameworks?

They are built on fundamentally different objects — CI is a wavefunction method expanding in explicit many-electron basis functions, while DFT works with the electron density and an (in principle exact, in practice approximate) exchange-correlation functional — so they are not directly combined in the same calculation. However, they are used together in practice: DFT orbitals are sometimes used as a starting point for building CI excitation spaces in multireference contexts, and full CI or high-level CISD/MRCI benchmarks are routinely used to validate and parametrize new density functionals against numerically exact reference data, especially for excited states and bond-breaking, where DFT's single-determinant-like starting point can struggle.

Before this

Where this leads