Integral-Direct Planewave Density Fitting

• Integral-direct planewave density fitting is a method that uses auxiliary plane-wave functions to represent and directly compute densities, turning costly nonlocal operations into efficient, diagonalizable tasks.
• It finds applications in electronic structure theory through planewave hybrid DFT and mixed Gaussian/plane-wave density fitting, which reduce computational bottlenecks and storage demands.
• In boundary-integral problems for Helmholtz and Maxwell equations, the technique regularizes singularities by replacing difficult kernel evaluations with on-the-fly, direct integral computations.

Searching arXiv for papers on integral-direct planewave density fitting and closely related methods. Integral-direct planewave density fitting denotes a family of methods in which densities, pair densities, or boundary densities are represented, interpolated, or contracted through plane-wave auxiliary functions while the governing operator is applied directly, rather than through stored four-index integrals or specialized singular-integration machinery. In electronic-structure theory, this idea appears in planewave hybrid density functional theory, mixed Gaussian/plane-wave density fitting for periodic Coulomb integrals, and all-electron periodic Fourier-transformed Coulomb methods (Hu et al., 2017, Sun et al., 2017, Dinh et al., 19 Sep 2025). In boundary-integral formulations for Helmholtz and Maxwell problems, closely related planewave density interpolation constructs local homogeneous solutions whose traces regularize weakly singular, strongly singular, hypersingular, and nearly singular kernels so that standard quadrature becomes applicable (Pérez-Arancibia et al., 2019, Pérez-Arancibia et al., 2019). Across these settings, the common structural feature is an auxiliary plane-wave representation that converts a difficult nonlocal operation into a diagonal or bounded reciprocal-space or quadrature problem.

1. Terminological scope and conceptual core

In the electronic-structure literature, “density fitting” or “resolution of identity” denotes an auxiliary-basis approximation to pair densities such as $$\rho_{\mu\nu}(\mathbf r)=\phi_\mu(\mathbf r)\phi_\nu(\mathbf r)$$ or to orbital products entering exchange and Coulomb operators. In the boundary-integral literature, “planewave density interpolation” or “planewave density fitting” refers to a local pointwise construction in which the density and, when needed, its normal derivative are matched by traces of plane-wave superpositions that solve the underlying PDE exactly (Pérez-Arancibia et al., 2019). The latter terminology emphasizes exact matching of prescribed tangential derivatives at a target point, whereas “fitting” can also include least-squares or pseudoinverse constructions for higher-order local expansions.

“Integral-direct” has a correspondingly specific meaning in each domain. In planewave hybrid DFT, standard exchange is already integral-direct in the sense that the action of the exchange operator is evaluated on the fly by FFT-based convolution rather than by storing four-center exchange integrals (Hu et al., 2017). In all-electron periodic Coulomb builds, integral-direct planewave density fitting means that plane-wave coefficients are obtained analytically by direct Fourier integrals of AO pair densities, so that no least-squares regression to a plane-wave basis and no storage of reciprocal-space density-fitting intermediates are required (Dinh et al., 19 Sep 2025). In boundary-integral equations, the method is integral-direct because the singular kernel is regularized at the integrand level and then evaluated with elementary quadrature, rather than by singularity extraction, Duffy transforms, or semi-analytic element integrals (Pérez-Arancibia et al., 2019).

A recurrent misconception is that these methods globally replace the primary discretization by plane waves. The boundary-integral formulations do not do this: the unknown surface current or density remains represented by RWG, Nyström, or Galerkin degrees of freedom, and the plane-wave construction is local and auxiliary (Pérez-Arancibia et al., 2019). Likewise, mixed Gaussian/plane-wave density fitting does not replace Gaussian AOs by plane waves; it partitions short-range and long-range structure so that compact Gaussian auxiliaries represent near-nuclear features while plane waves capture smooth periodic components (Sun et al., 2017).

2. Planewave hybrid DFT: integral-direct exchange, ISDF, and ACE

For hybrid DFT in a planewave basis, the Hartree–Fock-like exchange operator acting on a trial orbital $\psi_i$ can be written in real space as

$$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$

In a planewave code, this is evaluated integral-direct by forming the pair density $\rho_{ji}(r')=\phi_j(r')\psi_i(r')$, FFT-transforming it to reciprocal space, multiplying by the Coulomb kernel $v(G)=4\pi/G^2$, and inverse FFT-transforming to obtain the corresponding potential. For screened HSE06 exchange, the kernel is replaced by

$$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$

so the “Poisson solve” becomes a reciprocal-space convolution with $v_{\mathrm{SR}}(G)$ (Hu et al., 2017).

The computational bottleneck is that a conventional exchange application requires one Poisson-like or screened-Poisson solve per distinct pair density $\phi_j(r)\psi_i(r)$, giving $\mathcal O(N_e^2)$ FFT convolutions per inner diagonalization loop when $n\approx N_e$. Interpolative separable density fitting (ISDF) replaces this by the separable approximation

$\psi_i$0

where $\psi_i$1 are numerical auxiliary basis vectors and $\psi_i$2 are interpolation points selected by randomized QR with column pivoting on a compressed sample of the pair-product matrix. The interpolation-vector solve has the Galerkin form

$\psi_i$3

and the setup cost is $\psi_i$4 with a small preconstant. Because only $\psi_i$5 auxiliary functions are needed, with $\psi_i$6–$\psi_i$7 in practice, the number of Poisson or screened-Poisson convolutions is reduced from $\psi_i$8 to $\psi_i$9 (Hu et al., 2017).

The method is combined with the adaptively compressed exchange operator,

$$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$0

which is exact on the current trial-orbital subspace and reduces how often exchange must be applied. ISDF lowers the cost of each ACE update, while ACE reduces the number of updates to one per outer SCF iteration. The reported result is that ACE-ISDF reduces the computational cost associated with the exchange operator by nearly two orders of magnitude compared to existing approaches for a large silicon system with $$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$1 atoms; converged hybrid functional calculation results for a 1000-atom bulk silicon were obtained within 10 minutes on 2000 computational cores, and the method scaled to 8192 computational cores for a 4096-atom bulk silicon system (Hu et al., 2017).

The same study also reports that the method produces accurate energies and forces for insulating and metallic systems. Total energy error per atom and maximum force error can be kept below $$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$2 Hartree/atom and $$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$3 Hartree/Bohr with $$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$4–$$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$5, and energies, band gaps, and exchange energies converge systematically as $$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$6 increases. A plausible implication is that the method’s utility is not limited to wide-gap systems, because the required rank parameter is reported to be largely insensitive to the band gap (Hu et al., 2017).

3. Mixed Gaussian/plane-wave density fitting in periodic electronic structure

The mixed density fitting formulation for periodic systems uses both Gaussian and plane-wave auxiliary functions to approximate AO pair densities in the Coulomb metric (Sun et al., 2017). For a unit cell of volume $$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$7, the plane-wave auxiliary functions are

$$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$8

with reciprocal-space Coulomb kernel

$$(K\psi_i)(r)=-\sum_j^{\text{occ}}\phi_j(r)\int \frac{\phi_j(r')\psi_i(r')}{|r-r'|}\,dr'.$$9

The Fourier transform of the pair density is

$\rho_{ji}(r')=\phi_j(r')\psi_i(r')$0

so the three-index term and plane-wave metric become

$\rho_{ji}(r')=\phi_j(r')\psi_i(r')$1

This diagonal reciprocal-space metric is one of the principal attractions of a plane-wave auxiliary representation.

The Gaussian part of the mixed auxiliary basis is modified by smooth compensating Gaussians so that the periodic Gaussian fitting functions carry zero net charge and zero multipoles. This removes the problematic average-density component and allows the net $\rho_{ji}(r')=\phi_j(r')\psi_i(r')$2 term to be handled consistently under periodic boundary conditions. The Coulomb-metric orthogonalization then removes linear dependencies between Gaussian and plane-wave subspaces while preserving a diagonal plane-wave block. Eigenvectors below a threshold of $\rho_{ji}(r')=\phi_j(r')\psi_i(r')$3 are removed for numerical stability (Sun et al., 2017).

The mixed construction is motivated by the all-electron periodic setting. Near nuclei, steep Gaussian auxiliaries efficiently represent cusp-like short-range structure, whereas in interstitial regions the density is smooth and can be represented systematically by plane waves. The method therefore avoids the prohibitively high cutoffs that a pure plane-wave fitting basis would require for all-electron core structure. Benchmarks reported for a $\rho_{ji}(r')=\phi_j(r')\psi_i(r')$4-point hydrogen crystal and for all-electron silicon show that MDF can be substantially more accurate than pure plane-wave density fitting with the same number of plane waves; for the hydrogen crystal, MDF with 729 plane waves achieved approximately $\rho_{ji}(r')=\phi_j(r')\psi_i(r')$5 $\rho_{ji}(r')=\phi_j(r')\psi_i(r')$6 fitting error and was $\rho_{ji}(r')=\phi_j(r')\psi_i(r')$7–$\rho_{ji}(r')=\phi_j(r')\psi_i(r')$8 orders of magnitude more accurate than pure plane-wave DF with the same $\rho_{ji}(r')=\phi_j(r')\psi_i(r')$9 (Sun et al., 2017).

This framework is integral-direct in the usual density-fitting sense: ERIs are not stored as four-index objects but reconstructed on demand from lower-rank quantities. The plane-wave block is especially favorable because reciprocal-space contractions scale as $v(G)=4\pi/G^2$0 per contraction and exploit the diagonal Coulomb metric. The paper also emphasizes that analytical Fourier transforms, rather than discrete FFT sampling of steep Gaussians, are essential for all-electron accuracy; using discrete FFTs for steep Gaussians yields large errors and can destroy positive definiteness of the Coulomb metric (Sun et al., 2017).

4. All-electron periodic ID-PWDF and its relation to integral-direct periodic HF and MP2

The all-electron periodic Fourier-transformed Coulomb framework introduces an explicit “integral-direct planewave density fitting” long-range treatment inside an Ewald decomposition of the Coulomb matrix (Dinh et al., 19 Sep 2025). With $v(G)=4\pi/G^2$1,

$v(G)=4\pi/G^2$2

with the monopole $v(G)=4\pi/G^2$3 removed for charge neutrality. Short-range contributions are evaluated in real space using standard Gaussian density fitting, whereas long-range contributions are handled by plane waves.

The defining feature is that the plane-wave coefficients are not obtained by solving a least-squares problem in a reciprocal-space auxiliary basis. Instead, they are computed analytically by direct Fourier integrals of AO pair densities using the Gaussian product theorem and product-separable recurrences. For a primitive $v(G)=4\pi/G^2$4-Gaussian centered at $v(G)=4\pi/G^2$5 with exponent $v(G)=4\pi/G^2$6,

$v(G)=4\pi/G^2$7

For higher angular momentum, Hermite-polynomial recurrences are used. The long-range Coulomb contribution is then accumulated on the fly through

$v(G)=4\pi/G^2$8

$v(G)=4\pi/G^2$9

$$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$0

No plane-wave density-fitting intermediates are stored (Dinh et al., 19 Sep 2025).

The method further classifies contracted shells into compact and diffuse primitives using

$$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$1

Compact–compact and compact–diffuse terms use the Ewald short-range plus long-range path, while diffuse–diffuse terms are treated entirely by PW-DF in a GPW-style real-space density evaluation and FFT step. The reported timing data show orders-of-magnitude speedups over RSDF for dense $$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$2-point meshes. For the benzene crystal with cc-pVDZ, the QCPBC ID-PWDF RI-J time was 22.60 s at $$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$3, 28.61 s at $$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$4, 32.08 s at $$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$5, 33.11 s at $$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$6, and 42.18 s at $$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$7, whereas the closely related PySCF RSDF times were 22.86 s, 162.56 s, 840.61 s, 2153.08 s, and not reported at $$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$8 because PySCF was prohibitive (Dinh et al., 19 Sep 2025).

Accuracy is reported at the density-fitting level. For silicon with cc-pVDZ/TZ/QZ, SCF energy variation across $$v_{\mathrm{SR}}(r-r')=\frac{\operatorname{erfc}(\omega |r-r'|)}{|r-r'|}, \qquad v_{\mathrm{SR}}(G)=\frac{4\pi}{G^2}\Bigl(1-e^{-G^2/(4\omega^2)}\Bigr),$$9 is less than 20 $v_{\mathrm{SR}}(G)$0Ha up to $v_{\mathrm{SR}}(G)$1, and errors decrease with denser $v_{\mathrm{SR}}(G)$2-meshes while staying within typical DF errors of 50–60 $v_{\mathrm{SR}}(G)$3Ha/atom. The method was applied with dispersion-corrected PBE to compute benzene-crystal cohesive energies and CO adsorption on MgO(001), with results stated to be in good agreement with existing literature (Dinh et al., 19 Sep 2025).

A closely related periodic Gaussian DF-HF/MP2 implementation does not use planewaves in the published code, but it explicitly formulates the corresponding plane-wave auxiliary-basis adaptation (Bintrim et al., 2022). In that formulation, the auxiliary functions are

$v_{\mathrm{SR}}(G)$4

the Coulomb metric is diagonal,

$v_{\mathrm{SR}}(G)$5

and the three-index quantities are the Fourier components of the AO pair density,

$v_{\mathrm{SR}}(G)$6

The paper identifies this as the natural reciprocal-space analogue of integral-direct periodic DF, with on-the-fly blocking strategies for $v_{\mathrm{SR}}(G)$7, $v_{\mathrm{SR}}(G)$8, and MP2 contractions and with the $v_{\mathrm{SR}}(G)$9 singularity handled by setting the singular mode to zero together with an analytic correction such as a Madelung term (Bintrim et al., 2022). This suggests that the later all-electron ID-PWDF construction is part of a broader migration of periodic integral-direct DF toward analytically evaluated reciprocal-space auxiliaries.

5. Planewave density interpolation for Helmholtz and EFIE operators

In three-dimensional Helmholtz boundary integral equations, planewave density interpolation regularizes singular kernels by matching the density and its normal derivative with traces of local homogeneous Helmholtz solutions (Pérez-Arancibia et al., 2019). For the free-space Green’s function

$\phi_j(r)\psi_i(r)$0

the single- and double-layer potentials and the Calderón operators inherit weakly singular or hypersingular kernels. Green’s third identity is used in add–subtract form: for a homogeneous solution $\phi_j(r)\psi_i(r)$1, one inserts $\phi_j(r)\psi_i(r)$2 and its normal derivative against the true density so that the difference $\phi_j(r)\psi_i(r)$3 or $\phi_j(r)\psi_i(r)$4 cancels the singular behavior at $\phi_j(r)\psi_i(r)$5.

The local interpolant is expanded in plane waves,

$\phi_j(r)\psi_i(r)$6

and trace-matching conditions enforce, for all multi-indices $\phi_j(r)\psi_i(r)$7 with $\phi_j(r)\psi_i(r)$8,

$\phi_j(r)\psi_i(r)$9

Consequently,

$\mathcal O(N_e^2)$0

which implies bounded integrands for Brakhage–Werner when $\mathcal O(N_e^2)$1 and for Burton–Miller when $\mathcal O(N_e^2)$2 (Pérez-Arancibia et al., 2019).

The same principle is extended to the electric field integral equation on PEC surfaces (Pérez-Arancibia et al., 2019). There the scattered field is represented as

$\mathcal O(N_e^2)$3

with single-layer-like and gradient-of-single-layer terms involving the surface current $\mathcal O(N_e^2)$4 and its surface divergence. PWDI constructs a local plane-wave superposition

$\mathcal O(N_e^2)$5

such that $\mathcal O(N_e^2)$6 vanishes to prescribed tangential order while $\mathcal O(N_e^2)$7 matches the target scalar density or the divergence of the RWG basis function. The regularized kernels are then bounded both on-surface and off-surface, including nearly singular interactions between close or overlapping components (Pérez-Arancibia et al., 2019).

In both Helmholtz and EFIE settings, the numerical consequence is that standard quadrature rules become applicable irrespective of the original singularity. For smooth surfaces with Chebyshev-based Nyström discretization, spectral convergence is observed, with far-field errors converging at approximately third order for $\mathcal O(N_e^2)$8–$\mathcal O(N_e^2)$9 and approximately fifth order for $n\approx N_e$0–$n\approx N_e$1 as $n\approx N_e$2 increases. For Galerkin BEM with piecewise linear basis functions on triangles, observed convergence is second order, consistent with $n\approx N_e$3 BEM, and PWDI allows direct quadrature on simple and composite surfaces, including touching and intersecting configurations (Pérez-Arancibia et al., 2019). In the EFIE setting, the method permits non-conformal meshes across composite surfaces, and it does not alter EFIE conditioning; Calderón-type preconditioning or loop–star decompositions remain relevant (Pérez-Arancibia et al., 2019).

6. Shared algorithmic patterns, numerical issues, and limitations

Several design patterns recur across these otherwise distinct literatures. First, the auxiliary plane-wave representation is used because it diagonalizes or regularizes the dominant operator: the Coulomb kernel becomes a diagonal reciprocal-space multiplier in periodic electronic-structure methods, while in boundary integrals the plane-wave interpolant is an exact homogeneous solution that cancels kernel singularities pointwise (Sun et al., 2017, Pérez-Arancibia et al., 2019). Second, the methods are integral-direct because they avoid the central storage bottleneck of the naive formulation: four-index exchange or Coulomb integrals are not stored in the planewave hybrid and periodic Gaussian settings, and singular element integrals are not isolated and tabulated in the boundary-integral setting (Hu et al., 2017, Dinh et al., 19 Sep 2025, Pérez-Arancibia et al., 2019).

Charge neutrality and treatment of the singular reciprocal mode are also universal numerical issues. In periodic Coulomb problems, the $n\approx N_e$4 contribution must be removed or compensated. Mixed Gaussian/plane-wave DF subtracts the unit-cell average density and combines the resulting constant with the electron–nuclear and nuclear–nuclear $n\approx N_e$5 contributions (Sun et al., 2017). The all-electron ID-PWDF method removes $n\approx N_e$6 in both short-range DF metrics and long-range reciprocal sums (Dinh et al., 19 Sep 2025). The periodic HF/MP2 framework similarly identifies the $n\approx N_e$7 singularity as the key mode requiring analytic correction for reliable thermodynamic-limit convergence (Bintrim et al., 2022).

The principal limitations are formulation-specific. In planewave hybrid DFT, ISDF setup is $n\approx N_e$8 and the stored interpolation vectors and corresponding potentials scale as $n\approx N_e$9; the method is therefore most beneficial when exchange application dominates total cost (Hu et al., 2017). In mixed Gaussian/plane-wave DF, diffuse Gaussians and low-energy plane waves can become nearly linearly dependent, so Coulomb-metric orthogonalization and threshold tuning are essential (Sun et al., 2017). In the all-electron periodic FTC method, current short-range DF storage still scales between quadratic and cubic with supercell size before sparsity saturates, and the short-range DF solve is $\psi_i$00 per SCF if the inverse is not precomputed (Dinh et al., 19 Sep 2025). In Helmholtz and EFIE applications, very high curvature, edge singularities, extremely close interactions, or gradient evaluations near the surface may require higher interpolation order, more directions, or mesh refinement; BM and EFIE conditioning issues remain distinct from kernel regularization (Pérez-Arancibia et al., 2019, Pérez-Arancibia et al., 2019).

The literature therefore presents integral-direct planewave density fitting not as a single algorithm but as a methodological class. In one branch, it accelerates hybrid DFT and periodic Coulomb builds by compressing or analytically evaluating pair-density information in reciprocal space. In another, it regularizes singular boundary operators by replacing the density locally with traces of plane-wave solutions. What unifies these branches is the use of plane-wave auxiliary structure to convert an otherwise expensive or singular integral operation into a numerically tractable direct evaluation.