Supercell ISDF

---

This is an idea from the 2026 paper Ab Initio Many Body Quantum Embedding and Local Correlation in Crystalline Materials using Interpolative Separable Density Fitting (https://arxiv.org/abs/2601.16379) by Yang and colleagues. I want to provide some more details that are left out of the paper for brevity but that are important when implementing this factorization. I will change some of the nomenclature of the paper and reuse some other.

We want to factorize

For this, we introduce an evenly spaced grid , which we index with and calculate

by using the ISDF style expression for the center Coulomb matrix :

which we found from the Interpolative Decomposition of the Co-density

The ID interpolation points (indexed by ) can be computed using a pivoted Cholesky Decomposition of the grid metric matrix

The ID gives us matrices of the form

which represent an overtedermined system of equations, since we generally have equations for many equations with , which we need to solve times.

For overdetermined systems, we can run into problems as there might be no solution that satisfies all of our equations. But we can find a linear least squares solution by using

which gives us the normal equation

This we derive by means of calculating the gradient of with respect to and setting it to zero, which in our case gives us

which in matrix representation is

This boils down to just an expression for each column of the columns in . The pruned grid metric is sometimes also denoted as and the grid vs. pruned grid metric (sometimes called Right Hand Side - RHS vector) is sometimes denoted as . I find it clearer to use the index labels, as it shows the underlying math expression for both metrics are identical.

In theory we could just compute using standard solving procedures and then use FFT to calculate , but we also avoid this by realizing that both and project into the column space, which means we can reformulate the equation

and plugging this into our integral expression for

where we can pull out as it has no dependence on the full real space grid indexed by and . We now need to evaluate

which we can do by using the convolution theorem that states that

i.e. a convolution the real-domain becomes a multiplication in the Fourier transformed version of the space. Here

We can use this be realizing that if we write our ISDF integral representation in a separated form, such that one variable is integrated out

and transform it to one "function-like" equation per point

it is just a convolution and can therefore be solved by

where denotes a hadamard product between the matrices. This only works as long as our grid we apply the transformation along (i.e. ) is evenly spaced! Since this is the case, we can find the Fourier transformed Coulomb function as by expanding the standard Coulomb function as a Fourier series over the reciprocal lattice vectors

which we can plug into Poissons equation

Here we also expanded the Dirac delta. Now we can let the Laplace operator act on the expanded version and it pulls down two times , so we get

where is the reciprocal space lattice vectors at grid point .

Now we want to directly evaluate the double sided equation

by using the convolution theorem on the function

We solve it algorithmically by calculating

Now we can apply Perseval's theorem about the integral of Fourier transformed functions and , which states that

where an upper marks complex conjugation and denotes the volume of our unit cell. Applied to our functions this yields our final expression for

and we do so algorithmically by evaluating the inner product

Note that we can just take the real part here as there is only one Gamma point and our orbitals are real-valued.

With available, we can exploit our Cholesky vectors from the start to calculate by using intermediates and triangular solves