Density Fitting (RI)

Density fitting (DF), also known as the resolution of the identity (RI) approximation, is the most widely used way to reduce the cost of the two-electron repulsion integrals (ERIs). The general idea is to decompose the four-center integrals into a contraction between two-center integrals, built with a set of pre-fitted auxiliary basis functions, and three-center integrals, built with the auxiliary basis and the original basis functions:

\[ (\mu\nu|\lambda\sigma) \approx \sum_{PQ}(\mu\nu|P)[\mathbf{V}^{-1}]_{PQ}(Q|\lambda\sigma) = \sum_P B^P_{\mu\nu} B^P_{\lambda\sigma}, \]

where \(P, Q\) index the auxiliary basis functions (fitting vectors), \(V_{PQ} \equiv (P|Q) = \int \phi_P(\mathbf{r}_1)\frac{1}{|\mathbf{r}_1-\mathbf{r}_2|}\phi_Q(\mathbf{r}_2)\,\mathrm{d}\mathbf{r}_1\mathrm{d}\mathbf{r}_2\) is the two-center integral formed from the auxiliary basis functions (the Coulomb metric), and \((Q|\lambda\sigma)\) is the three-center integral. The half-Coulomb-contracted three-center tensor \(B^P_{\mu\nu} = \sum_Q [\mathbf{V}^{-1/2}]_{PQ}(Q|\mu\nu)\) is formed once and reused. Physically, DF exploits the fact that products of basis functions (generalized charge densities) are nearly linearly dependent and can be represented compactly in an auxiliary basis whose size grows only linearly with the system.

The Coulomb metric is the standard and most robust choice, but other metrics (overlap, anti-Coulomb, attenuated Coulomb) are also possible and trade off accuracy against the locality/sparsity of the three-center tensor.

Scaling and when to use it

Density fitting does not change the asymptotic scaling of the calculation — generating the three-center integrals scales as \(O(N_{\mathrm{AO}}^2 N_{\mathrm{aux}})\), the Coulomb (J) build similarly, and the exchange (K) build as \(O(N_{\mathrm{AO}}^2 N_{\mathrm{occ}} N_{\mathrm{aux}})\). What it reduces is the prefactor: it avoids redundant four-center integral evaluations and replaces them with dense matrix multiplications (precompute-and-reuse of the lower-order tensors), which map extremely well onto GPU linear-algebra kernels.

As a result, DF is most advantageous for small-to-medium systems with large AO basis sets — roughly up to ~100 atoms or ~2000 basis functions. For sufficiently large systems the integral-direct exact-ERI approach, which can exploit asymptotic sparsity, eventually wins, because the dense linear algebra of DF does not benefit from that sparsity. The accuracy cost is negligible: in TeraChem, DF excitation energies agree with the integral-direct reference to within ~10 microhartree.

DF is especially impactful for CASSCF, where each macroiteration requires multiple J builds; TeraChem's GPU-accelerated DF-CASSCF (with full analytic gradients) accelerates total energy and gradient evaluations by roughly 3–30× for small-to-medium systems with large basis sets, making high-throughput multireference dynamics far more affordable.

Implementation

DF is implemented in both the J and K Fock-build routines and in the corresponding analytic-gradient routines. The two- and three-center integrals (and their nuclear derivatives) are computed with the McMurchie–Davidson algorithm in CUDA kernels generated by symbolic-mathematics meta-programming, with Schwarz pair screening (threshold \(10^{-11}\)) reducing the three-center work. For the exchange term, TeraChem provides two strategies: a density-matrix formulation (compatible with the standard J/K routines) and an MO-coefficient formulation that exploits the low rank of the density for better performance (with an SVD-based rank reduction for the generalized, non-positive-semidefinite densities that arise in CASSCF gradients and nonadiabatic couplings).

With the keyword use_density_fitting specified, TeraChem will automatically perform density fitting for the ERIs' Fock build and gradient routines.

The df_auxiliary_basis takes the same format as the basis keyword, either a TeraChem built-in basis set name or a path to a basis set file. Leaving it empty or setting it to default will use the default auxiliary basis set based on the basis keyword, provided the latter has a corresponding built-in auxiliary basis set.

Keywords

Keyword	Type	Default	Description
use_density_fitting	string	no	Turn on density fitting approximation for J and/or K routines (fock build and gradient). Valid options: j, k, jk/kj/yes, no
df_auxiliary_basis	string	default	Default, or the name of or the path to the auxiliary basis set to use for density fitting
df_metric_thre	float	1e-7	Threshold for density fitting metric matrix
df_decompress_pij	bool	no	Whether to decompress the pair-screened raw 3-center integrals
df_num_check	bool	no	Whether to perform numerical checks of the density fitting gradients
df_orb_coeff_k	bool	yes	Whether to use orbital coefficients for K
df_custom_kernel	bool	no	Whether to use custom kernel for density fitting K build (less performant)
df_int_direct_j	bool	no	Whether to use direct integrals for J
df_aux_true_spherical	bool	no	Whether to use spherical auxiliary basis functions

System-specific guidelines

The current implementation of density fitting stores the 3-center integrals, \(\widetilde{(P|\mu\nu)} = (P|Q)^{-1/2}(Q|\mu\nu)\), in the device memory. When the system is large, it is recommended to set df_orb_coeff_k to yes (the K routines will use the MO coefficients as input, resulting in a smaller intermediate buffer size) and df_decompress_pij to no (the pair-screened raw 3-center integrals has a smaller memory size).

Auxiliary basis set defaults

The following tables list the default auxiliary basis sets that will be used if the df_auxiliary_basis keyword is set to default.

Pople basis sets

Basis set	J-fit	JK-fit	RI-fit	OptRI
6-311G(d,p)	-	-	6-311gss-rifit	-
6-311G**	-	-	6-311gss-rifit	-
6-31G(d,p)	-	-	6-31gss-rifit	-
6-31G**	-	-	6-31gss-rifit	-

Dunning correlation-consistent basis sets

Basis set	J-fit	JK-fit	RI-fit	OptRI
cc-pVDZ	-	-	cc-pvdz-rifit	-
cc-pVTZ	-	cc-pvtz-jkfit	cc-pvtz-rifit	-
cc-pVQZ	-	cc-pvqz-jkfit	cc-pvqz-rifit	-
cc-pV5Z	-	cc-pv5z-jkfit	cc-pv5z-rifit	-
aug-cc-pVDZ	-	-	aug-cc-pvdz-rifit	aug-cc-pvdz-optri
aug-cc-pVTZ	-	cc-pvtz-jkfit	aug-cc-pvtz-rifit	aug-cc-pvtz-optri
aug-cc-pVQZ	-	cc-pvqz-jkfit	aug-cc-pvqz-rifit	aug-cc-pvqz-optri
aug-cc-pV5Z	-	cc-pv5z-jkfit	aug-cc-pv5z-rifit	aug-cc-pv5z-optri

Karlsruhe basis sets

Basis set	J-fit	JK-fit	RI-fit	OptRI
def2-SV(P)	def2-universal-jfit	def2-sv(p)-jkfit	def2-sv(p)-rifit	-
def2-SVP	def2-universal-jfit	def2-universal-jkfit	def2-svp-rifit	-
def2-TZVP	def2-universal-jfit	def2-universal-jkfit	def2-tzvp-rifit	-
def2-TZVPP	def2-universal-jfit	def2-universal-jkfit	def2-tzvpp-rifit	-
def2-QZVP	def2-universal-jfit	def2-universal-jkfit	def2-qzvp-rifit	-
def2-QZVPP	def2-universal-jfit	def2-universal-jkfit	def2-qzvpp-rifit	-

N.B., A dash (-) indicates that no auxiliary basis is available for that fitting type.

Reference

The GPU-accelerated density-fitting implementation in TeraChem and its application to CASSCF are described in:

R. Wang, Y. Wang, L. Lu, D. Hait, and T. J. Martínez, "Complete Active Space Self-Consistent Field with GPU-Accelerated Density Fitting," J. Chem. Theory Comput. 22, 3398–3414 (2026). doi:10.1021/acs.jctc.5c02079