Effective core potential theory

Effective core potential theory#

Heavy elements make explicit all-electron calculations expensive in two independent ways: there are simply more electrons to correlate/converge, and — for the heaviest elements — relativistic effects on the core electrons become significant even though chemistry happens almost entirely in the valence shell. An effective core potential (ECP), also called a pseudopotential, replaces the chemically inert core electrons with a fixed, parametrized one-electron operator, cutting both the electron count and the basis size while (for a well-fitted ECP) implicitly folding in scalar relativistic effects the core electrons would otherwise need explicit treatment for. This chapter derives the semilocal ECP operator qc-rs implements, its non-obvious radial-power convention, and the semi-numerical evaluation method that makes it numerically robust for heavy, tightly contracted primitives — grounded in .design/24qc.scf-ecp.md and qc-rs’s AGENTS.md qc-ecp rules.

The valence-only Hamiltonian#

With nuclear charge \(Z_A\) and \(N_A^{\text{core}}\) core electrons removed by the ECP, the effective nuclear charge is \(Z_A^{\text{eff}}=Z_A-N_A^{\text{core}}\), and the non-relativistic valence-only electronic Hamiltonian becomes

\[ \hat H = \sum_i\left[-\tfrac12\nabla_i^2 - \sum_A\frac{Z_A^{\text{eff}}}{r_{iA}} + \sum_A\hat U_A^{\text{ECP}}(i)\right] + \sum_{i<j}\frac{1}{r_{ij}} + E_{NN}^{\text{eff}}, \qquad E_{NN}^{\text{eff}} = \sum_{A<B}\frac{Z_A^{\text{eff}}Z_B^{\text{eff}}}{R_{AB}}. \]

Every piece of this — the reduced nuclear attraction, the reduced nuclear repulsion, and the new \(\hat U_A^{\text{ECP}}\) operator itself — flows from a single substitution: an all-electron atom is simply the special case \(N_A^{\text{core}}=0\), \(\hat U_A^{\text{ECP}}=0\). The valence electron count follows directly, \(N_e=\sum_AZ_A^{\text{eff}}-Q\) (with \(Q\) the molecular charge), and this is the electron count qc-rs’s spin/multiplicity parity checks use — an ECP calculation’s electron-counting logic is not a special case bolted onto the all-electron logic, it is the all-electron logic applied to \(Z_A^{\text{eff}}\) instead of \(Z_A\).

The semilocal scalar ECP operator#

The chemically important part is \(\hat U_A^{\text{ECP}}\) itself. A semilocal ECP (the scalar/AREP — averaged relativistic effective potential — variety qc-rs implements) splits into a purely radial local channel plus a sum of angular-momentum projectors, each carrying its own radial correction:

\[ \hat U_A^{\text{scalar}} = U_{A,L}^{\text{loc}}(r_A) + \sum_{l=0}^{L-1}\sum_{m=-l}^{l} |Y_{lm}^A\rangle\,\Delta U_{A,l}(r_A)\,\langle Y_{lm}^A|, \qquad \Delta U_{A,l}(r) = U_{A,l}(r) - U_{A,L}^{\text{loc}}(r), \]

where \(Y_{lm}^A\) are real spherical harmonics centered on atom \(A\) and \(L\) is the highest angular momentum explicitly represented in the ECP (electrons with \(l\ge L\) see only the local channel — the physical picture being that beyond some angular momentum, the core no longer needs a channel-specific correction). The projector \(|Y_{lm}^A\rangle\langle Y_{lm}^A|\) means an electron of angular momentum \(l\) relative to center \(A\) experiences \(U_{A,L}^{\text{loc}}+\Delta U_{A,l}\), while any other electron experiences only the local channel — this is exactly the physical statement that the pseudopotential’s channel-dependence represents the different degree to which the Pauli exclusion principle (orthogonality to the removed core orbitals) affects electrons of different angular momentum. Because published ECP data sometimes tabulates the local channel directly as ul rather than requiring the consumer to compute \(\Delta U\) from raw per-channel values, qc-rs does not silently convert to a difference potential on input — it stores exactly what the database provides and applies the subtraction (if needed) only at evaluation time.

Each radial channel is itself a sum of primitive Gaussian-type radial functions,

\[ U_{A,l}(r) = \sum_k C_{A,lk}\,r^{n_{A,lk}-2}\,e^{-\alpha_{A,lk}r^2}. \]

The radial-power convention is a documented trap, not a cosmetic detail. The stored integer \(n_{A,lk}\) — called RADI_POWER in NWChem/PySCF/libecpint file conventions and EcpPrim::n in qc-mol — is kept and passed through exactly as read from the basis file, with no -2 correction applied in qc-rs’s own code, even though the formula above visibly has an \(r^{n-2}\) dependence. The -2 is already baked into the file convention itself (verified against PySCF’s ECPscalar); silently “fixing” it by subtracting 2 again inside qc-rs would double-apply the correction and produce a wrong radial power. This is precisely the kind of convention mismatch the basis-and-AO chapter warns about more generally: a stored integer’s meaning is defined by the file format that produced it, not by what a textbook formula’s exponent happens to look like.

The full AO matrix element is a straightforward one-electron expectation value over every ECP center:

\[ V^{\text{ECP}}_{\mu\nu} = \sum_A\langle\chi_\mu|\hat U_A^{\text{ECP}}|\chi_\nu\rangle. \]

Since \(U_A^{\text{ECP}}\)’s projector structure is a genuinely semilocal operator (it is not simply \(1/r\) or a Gaussian multiplying the identity — it depends on the angular momentum of the electron relative to center \(A\), which is not diagonal in an arbitrary AO basis centered elsewhere), it falls entirely outside libcint’s operator algebra, which is why qc-ecp needs its own kernel rather than reusing the production libcint one/two-electron integral drivers the rest of qc-rs’s Hamiltonian assembly is built on.

Primitive normalization#

The ECP radial primitive’s normalization, verified against PySCF’s ECPscalar, is

\[ N(l,\zeta) = \sqrt{\frac{2\,(2\zeta)^{l+3/2}}{\Gamma(l+3/2)}}, \]

matching a unit radial self-overlap for the primitive — structurally the same normalization idea as ordinary Gaussian AO primitives (derived in the basis-and-AO chapter), just specialized to the ECP radial-power convention above rather than the ordinary GTO angular-momentum exponent. Because the full AO matrix reuses qc-grid’s GaussianShell and spherical-value evaluator (rather than a bespoke ECP-only AO representation), normalization, general contraction, and AO ordering all match libcint’s own conventions automatically — the ECP matrix is built directly in libcint AO order, with no separate reordering pass needed downstream.

Evaluation: analytic angular factors, semi-numerical radial quadrature#

The semilocal projector’s angular dependence has a genuine closed-form analytic treatment — real spherical-harmonic angular factors, Gaunt-coefficient-factorized (the same Gaunt/Clebsch-Gordan machinery that couples angular momenta when multiplying spherical harmonics) — but the radial integral over the product of two (possibly widely separated, off-center) Gaussian AO primitives and the ECP’s own radial Gaussian expansion has no similarly clean closed form for the general off-center case. qc-rs’s production path is therefore explicitly semi-numerical: exact angular factors, numerical radial quadrature — the same general strategy PySCF/libecpint use, not a qc-rs-specific compromise.

The production driver, qc-ecp::grid_radial::ecp_matrix_grid, is a direct algorithmic port of PySCF’s nr_ecp.c (independently reimplemented in Rust, reproducing its numbers): a nested Gauss-Chebyshev radial grid with per-shell-pair adaptive refinement. Rather than fixing one radial point count for every shell pair regardless of how sharply peaked the integrand is, the grid starts coarse (31 points, “LEVEL0”) and doubles for a given shell-pair/channel block only until that block’s radial integral has converged to a tight tolerance (\(10^{-12}\)), up to a master grid of 2047 points (“LEVEL_MAX”):

Algorithm 6 (Adaptive radial-refinement ECP matrix element)

Input: a shell-pair/channel block’s radial integrand (a product of Gaussian radial functions and the ECP’s own radial expansion), convergence tolerance \(\tau=10^{-12}\). Output: the converged radial contribution to that block of the ECP matrix.

  1. Start at the coarsest refinement level \(L_0\) (31 master-grid points).

  2. Evaluate the radial estimate \(E_L = w_{\text{scale},L}\sum_{i\in\text{subset}_L}w_i\,f(r_i)\), a strided subsample of a single master Gauss-Chebyshev grid on \((0,\infty)\), with a level-dependent rescaling \(w_{\text{scale},L}=2^{L_{\max}-L}\) so every level estimates the same integral.

  3. If \(|E_L-E_{L-1}|<\tau\) (or \(L=L_0\), no prior estimate to compare), accept \(E_L\).

  4. Otherwise double the point count (advance to \(L+1\), a denser strided subset of the same master grid — no new grid is built) and repeat.

This adaptivity is what keeps the method both robust and fast simultaneously: smooth, well-separated shell pairs (chemically the common case — most AO pairs in a large molecule involve at least one function far from a given ECP center) converge at the coarsest level and cost only \(\sim31\) radial points, while only the sharp, tightly-contracted, close-together pairs — the ones that would actually need high radial resolution — pay for the expensive 2047-point refinement. A fixed-grid approach would have to use the expensive grid everywhere to be safe for the hard cases, or risk silent inaccuracy on them; adaptive refinement gets the accuracy of the dense grid only where the physics actually demands it.

Why this exists instead of, and alongside, a closed-form analytic radial path. qc-rs also has an analytic off-center radial path (matrix::ecp_matrix_analytic, via modified spherical Bessel functions and the same Gaunt factorization) that is exact in principle but has a documented numerical failure mode for tight heavy-element primitives: the closed-form expression splits into a prefactor that can itself diverge (\(e^{k^2/4p}\to\infty\) for large exponents) multiplied by a term that vanishes (\(e^{-\zeta d^2}\to0\) for a large separation-dependent exponent), and floating-point arithmetic evaluates \(\infty\times0\) as NaN rather than recovering the finite product the two factors’ combination should give analytically. The grid-based radial form sidesteps this because its integrand, \(e^{-\zeta(r-d)^2}\cdot M_\lambda(2\zeta dr)\) (a modified spherical Bessel function argument built directly from the combined, already-finite exponent rather than from two separately-exploding/vanishing factors), stays numerically bounded throughout — this is why ecp_matrix_grid rather than the analytic path is qc-rs’s production default for the general case, even though the analytic path remains available (and useful) for cross-validation.

The angular part carries over unchanged from the analytic path regardless of which radial method is used — it is validated once, independent of radial strategy, and reused: for the “Type 1” (single-center, same-atom AO pair) and “Type 2” (the AO pair’s centers differ from each other, or from the ECP center, or both) angular structures, both a Cartesian-monomial-based path and the spherical-harmonic/Gaunt path are implemented and cross-validated against each other to \(\sim10^{-12}\)–\(10^{-13}\), with the faster one selected as the default for each case (Cartesian for Type 1 at high angular momentum, Gaunt for Type 2 where its per-center caching wins).

Verified example: a heavy element with an ECP basis#

import qc
m = qc.chk.new(atom="I 0 0 0; I 0 0 2.666", ao="lanl2dz", unit="angstrom").scf(ref="r").run()
m.scf.energy, m.scf.converged
# (-22.307234220033017, True)
m.scf.energy_components
# {'core': -56.253929996038075, 'coulomb': 28.188845752921235, 'exchange': -3.968212743262761}

I₂ has \(Z=53\) per atom (106 electrons all-electron), but LANL2DZ’s ECP removes the bulk of each iodine’s core, leaving only the valence electrons (\(Z^{\text{eff}}\) per atom, far fewer than 53) to appear explicitly in the SCF. The resulting total electronic energy, \(\sim-22.3\ E_h\), is far smaller in magnitude than a genuine all-electron iodine dimer’s would be — exactly the expected signature of an ECP calculation: the deep, chemically inert core energy has been folded into the fixed \(\hat U_A^{\text{ECP}}\) operator and never appears as an explicit energy contribution at all. The core component above already includes both the reduced nuclear attraction and the ECP operator’s own contribution — qc-rs does not report a separate “ECP energy” term, because \(\hat U_A^{\text{ECP}}\) is architecturally just another piece of the valence-only core Hamiltonian \(T+V_{ne}[Z^{\text{eff}}]+\hat U^{\text{ECP}}\), not a bolt-on correction computed and reported separately.

Exercise 13

  1. The radial-power convention keeps RADI_POWER exactly as read from the basis file, with no -2 correction applied in qc-rs’s own code, even though the radial formula visibly contains \(r^{n-2}\). Explain in one sentence why “helpfully” subtracting 2 again inside qc-rs, to make the code’s exponent match the textbook formula literally, would actually be a bug rather than a clarification.

  2. The adaptive radial grid doubles the point count only for shell-pair/channel blocks that have not yet converged, rather than using the finest (2047-point) grid uniformly. What specific class of shell pairs would make a uniform coarse grid dangerous, and what would make a uniform fine grid wasteful?

  3. The closed-form analytic ECP radial integral can literally evaluate to NaN for tight heavy-element primitives, even though the true mathematical value is finite. Identify the two colliding limits (one diverging, one vanishing) responsible, and explain why reorganizing the same physics into the grid-based integrand avoids the problem rather than merely working around it with clamping or a special case.

The next chapter, on implicit solvation, covers a second Hamiltonian modification qc-rs supports — but one that, unlike the ECP’s fixed operator, is genuinely density-dependent and must therefore participate in the SCF self-consistency loop and, per the linear-response chapter, in the response kernel as well.