Two-electron integrals can be decomposed using Cholesky
decomposition^{Krylov:2013b} giving rise to the same representation as in
RI and substantially reducing the cost of integral transformation, disk storage
requirements, and improving parallel performance:

$$(\mu \nu |\lambda \sigma )\approx \sum _{P=1}^{M}{B}_{\mu \nu}^{P}{B}_{\lambda \sigma}^{P},$$ | (6.39) |

The rank of Cholesky decomposition, $M$, is typically 3-10 times larger than
the number of basis functions $N$ (Ref. Aquilante:2009); it
depends on the decomposition threshold $\delta $ and is considerably smaller
than the full rank of the matrix,
$N(N+1)/2$ (Refs. Aquilante:2009, Beebe:1977, Wilson:1990).
Cholesky decomposition removes linear dependencies in product
densities $(\mu \nu |$,^{Aquilante:2009} allowing one to obtain compact
approximation to the original matrix with accuracy, in principle, up to machine
precision.

Decomposition threshold $\delta $ is the only parameter that controls accuracy and the rank of the decomposition. Cholesky decomposition is invoked by specifying CHOLESKY_TOL that defines the accuracy with which decomposition should be performed. For most calculations tolerance of $\delta ={10}^{-3}$ gives a good balance between accuracy and compactness of the rank. Tolerance of $\delta ={10}^{-2}$ can be used for exploratory calculations and $\delta ={10}^{-4}$ for high-accuracy calculations. Similar to RI, Cholesky-decomposed integrals can be transformed back, into the canonical MO form, using CC_DIRECT_RI keyword.