Exact 4-center ERIs

The default way TeraChem handles the two-electron repulsion integrals (ERIs) is to evaluate the exact four-center integrals directly on the GPU, recomputing them on the fly every SCF iteration rather than storing them. This integral-direct strategy is what made TeraChem's original GPU acceleration possible, and it remains the most accurate and robust option. The density-fitting and THC pages describe approximations to this baseline.

Integral-direct SCF

For a system with more than a few hundred basis functions, storing the full ERI tensor is infeasible (it would require tens of terabytes for ~1000 basis functions), and transferring integrals between GPU and CPU is slower than recomputing them. TeraChem therefore builds the Coulomb and exchange matrices entirely on the GPU from freshly evaluated primitive integrals, and only the resulting matrices (size \(O(N^2)\)) are returned to the CPU. The Fock matrix is then

\[ \mathbf{F}(\mathbf{C}) = \mathbf{H}_{\mathrm{core}} + \mathbf{J}(\mathbf{C}) - \tfrac{1}{2}\mathbf{K}(\mathbf{C}) , \]

with

\[ J_{\mu\nu} = \sum_{\lambda\sigma}(\mu\nu|\lambda\sigma)D_{\lambda\sigma} , \qquad K_{\mu\nu} = \sum_{\lambda\sigma}(\mu\lambda|\nu\sigma)D_{\lambda\sigma} . \]

The J and K builds use different GPU algorithms because they have different memory-access patterns: the J build is formed directly from primitive integrals in a Hermite-Gaussian (McMurchie–Davidson) representation with the density matrix pre-contracted, while the K build assigns one thread block per matrix element. Both are mapped so that each GPU thread evaluates one primitive integral (or a small batch for higher angular momentum), with separate hand- and machine-generated kernels for each angular-momentum class (ss|ss, ss|sp, sp|sp, ss|pp, sp|pp, pp|pp, …). With efficient screening the Fock build scales empirically close to \(O(N^2)\), even though the formal count of integrals is \(O(N^4)\).

Schwarz screening

TeraChem prescreens integrals with the Schwarz inequality, which bounds each ERI by the product of the bra and ket pair quantities,

\[ |(\mu\nu|\lambda\sigma)| \le \sqrt{(\mu\nu|\mu\nu)}\,\sqrt{(\lambda\sigma|\lambda\sigma)} . \]

The bra and ket pairs are sorted by their Schwarz upper bound so that negligible integrals are grouped together and skipped, which is what reduces the effective scaling. The threshold below which an integral is discarded is set by threall (gradient contributions use thregr). In the K build, the bound is additionally weighted by the density-matrix elements (with a small guard parameter) so that density sparsity is exploited as well.

Numerical precision

Because GPUs deliver much higher throughput in single precision, precision is an important control. TeraChem evaluates ERIs in single, double, or mixed precision, accumulating the largest contributions in double precision to control error. The behavior is selected with the precision keyword:

`precision`	Meaning
`single`	All integrals in 32-bit. Fastest, least accurate; appropriate mainly for the early SCF iterations or rough work.
`mixed`	Integrals above the `threspdp` magnitude (default \(10^{-5}\)) are computed in double precision and the rest in single. A good speed/accuracy compromise and the usual recommendation.
`double`	All integrals in 64-bit. Most accurate and slowest; needed for tight gradients, frequencies, and high-accuracy energetics.
`dynamic`	Like `mixed`, but the single/double cutoff is tightened automatically as the SCF converges. (Not supported for CASSCF.)

Even when the integrals themselves are computed in single precision, the J and K matrix elements are accumulated in double precision to avoid error build-up from summing many integrals of widely varying magnitude. An incremental Fock build (forming only the change in J/K from the density difference between successive iterations) further improves accuracy with limited-precision integrals and reduces work as the SCF converges.

Keywords

Keyword	Type	Default	Description
`precision`	string	`mixed` (typical)	ERI precision strategy: `single`, `mixed`, `double`, or `dynamic` (see table above).
`threspdp`	float	\(10^{-5}\) (mixed)	Magnitude threshold separating single- from double-precision integrals in `mixed`/`dynamic` precision.
`threall`	float	internal	Schwarz screening threshold for ERIs in the Fock build; integrals bounded below this are skipped.
`thregr`	float	internal	Screening threshold for the gradient (derivative-integral) contributions.
`xtol`	float	internal	Tolerance controlling basis-pair/grid screening.

Gradients

Analytic nuclear gradients require the derivative ERIs, which are evaluated with the same GPU integral engine (generalized derivative J and K matrices) and benefit from the same screening and precision controls. Tight gradients and vibrational frequencies generally warrant precision double.

References

The GPU integral-direct engine is described in the "Quantum Chemistry on Graphical Processing Units" series:

I. S. Ufimtsev and T. J. Martínez, "Quantum Chemistry on Graphical Processing Units. 1. Strategies for Two-Electron Integral Evaluation," J. Chem. Theory Comput. 4, 222–231 (2008). doi:10.1021/ct700268q
I. S. Ufimtsev and T. J. Martínez, "Quantum Chemistry on Graphical Processing Units. 2. Direct Self-Consistent-Field Implementation," J. Chem. Theory Comput. 5, 1004–1015 (2009). doi:10.1021/ct800526s
I. S. Ufimtsev and T. J. Martínez, "Quantum Chemistry on Graphical Processing Units. 3. Analytical Energy Gradients, Geometry Optimization, and First Principles Molecular Dynamics," J. Chem. Theory Comput. 5, 2619–2628 (2009). doi:10.1021/ct9003004