Non-Uniform Fast Fourier Transform for Two-Electron Repulsion Integrals

Two-electron repulsion integrals (ERIs) constitute a fundamental computational bottleneck in quantum chemistry, scaling as $O(N^4)$ with the number of basis functions. This Projects presents implementation and analysis of ERI computation using the Non-Uniform Fast Fourier Transform (NUFFT) on adaptive Becke grids. We develop a modular six-component framework encompassing Gaussian-type orbital evaluation, numerical grid generation, orbital product computation, Fourier transformation, two-electron integral calculation, and systematic method comparison. Numerical validation demonstrates that all basis functions (s, p, d types) satisfy normalization to within $10^{-15}$ to $10^{-10}$ relative error, overlap integrals achieve $3.16 \times 10^{-10}$ accuracy against analytical values, and NUFFT-based Fourier transforms agree with analytical results to $10^{-12}$. For ERI computation, NUFFT on Becke grids achieves a relative error of $4.54 \times 10^{-6}$ using 7,760 grid points, while conventional FFT on uniform grids yields $1.05 \times 10^{-1}$ relative error (10.5%) with 884,736 points. This corresponds to a 23,038-fold improvement in accuracy using 114 times fewer grid points. Direct real-space integration fails with errors exceeding $10^{6}$ due to the Coulomb singularity at $r_{12} = 0$, which Fourier methods handle analytically through the $4\pi/k^2$ kernel. All results are validated against analytical reference values computed via the Boys function.

NUFFT achieves 23,038× better accuracy with 114× fewer grid points compared to FFT.

Two-electron repulsion integrals are the computational bottleneck of quantum chemistry, scaling as $O(N^4)$:

$$(\mu\nu|\sigma\lambda) = \int!!\int \frac{\chi_\mu(\mathbf{r}_1)\chi_\nu(\mathbf{r}_1)\chi_\sigma(\mathbf{r}_2)\chi_\lambda(\mathbf{r}_2)}{|\mathbf{r}_1 - \mathbf{r}_2|} , d\mathbf{r}_1 , d\mathbf{r}_2$$

The Coulomb operator $1/|\mathbf{r}_1 - \mathbf{r}2|$ has a singularity at $r{12}=0$ that destroys direct numerical integration. The key insight is the Fourier representation:

$$\frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} = \frac{1}{2\pi^2} \int \frac{e^{i\mathbf{k}\cdot(\mathbf{r}_1 - \mathbf{r}_2)}}{k^2} , d\mathbf{k}$$

which converts the 6D real-space integral into a 3D Fourier-space integral with the well-behaved $4\pi/k^2$ kernel:

$$(\mu\nu|\sigma\lambda) = \frac{1}{2\pi^2} \int \frac{\tilde{\rho}_{\mu\nu}(\mathbf{k}) \cdot \tilde{\rho}^*_{\sigma\lambda}(\mathbf{k})}{k^2} , d\mathbf{k}$$

Implements Cartesian Gaussian-Type Orbitals for s-type ($l=0$), p-type ($l=1$), and d-type ($l=2$) angular momenta, including contracted basis functions (STO-3G).

A Cartesian GTO centered at $\mathbf{R}$ is:

$$\chi(\mathbf{r}) = N \cdot (x - R_x)^{i_x} (y - R_y)^{i_y} (z - R_z)^{i_z} \cdot e^{-\alpha|\mathbf{r}-\mathbf{R}|^2}$$

with normalization constant:

$$N = \left(\frac{2\alpha}{\pi}\right)^{3/4} \cdot (4\alpha)^{l/2} \cdot \frac{1}{\sqrt{(2i_x-1)!! \cdot (2i_y-1)!! \cdot (2i_z-1)!!}}$$

Primitive GTO evaluation — vectorized over all grid points simultaneously:

1. Compute displacement vectors $\mathbf{d} = \text{points} - \text{center}$
2. Compute $r^2 = d_x^2 + d_y^2 + d_z^2$
3. Radial part: exp(-α * r²)
4. Angular part: $d_x^{i_x} \cdot d_y^{i_y} \cdot d_z^{i_z}$
5. Return $N \cdot \text{angular} \cdot \text{radial}$

Contracted GTO (STO-3G) — sums $K$ primitives: $\chi = \sum_{k=1}^K c_k N_k e^{-\alpha_k r^2}$

Shell evaluation — allocates output array of shape (ncomp, npts) where $n_\text{comp} = (l+1)(l+2)/2$, fills each row with the contracted GTO for the corresponding $(i_x, i_y, i_z)$ tuple (enumerated by looping $i_x$ from $l$ down to 0, $i_y$ from $l - i_x$ down to 0, $i_z = l - i_x - i_y$).

Double factorial — iterative: start at $k = 2n-1$, multiply into accumulator, decrement by 2 until $k \leq 0$.

Test 1.0 — Double Factorial Verification

Test 1.1 — GTO Normalization (Grid: $50^3 = 125{,}000$ pts, $L=20$ Bohr, $h=0.4082$ Bohr)

Test 1.2 — Gaussian Decay ($\alpha=1.0$, $N=0.712705$)

Confirms $e^{-\alpha r^2}$ behavior; negligible beyond 4 Bohr.

Test 1.3 — Angular Nodal Structure

Test 1.4 — Contracted GTO (STO-3G Hydrogen)

$\int|\chi_\text{STO-3G}|^2 d\mathbf{r} = 0.999984$ (finite grid error). Maximum at origin.

Test 1.5 — d-Shell Orthogonality

The overlap of $1/3$ between e.g. $d_{x^2}$ and $d_{y^2}$ is expected — Cartesian $x^2 e^{-\alpha r^2}$ and $y^2 e^{-\alpha r^2}$ are not the orthogonal spherical harmonic d-orbitals. The value arises from $\int x^2 y^2 e^{-2\alpha r^2} d\mathbf{r} \neq 0$.

Constructs two types of quadrature grids: uniform Cartesian grids (FFT-compatible) and adaptive atom-centered Becke grids (for NUFFT).

Uniform grid — fills a cubic box $[-L/2, L/2]^3$ with spacing $h = L/N$:

• 1D coordinates via linspace(-L/2, L/2, N, endpoint=False) (endpoint excluded for FFT periodicity)
• 3D tensor product via meshgrid → $N^3$ total points, all weights = $h^3$
• Constructor variants: from_box, from_center, from_molecule (auto-sizes around atomic coordinates with padding)

Becke grid — atom-centered spherical grid:

1. Radial: Gauss-Legendre points ${x_i, w_i^{GL}}$ on $[-1,1]$, mapped via Becke transform $r_i = R(1+x_i)/(1-x_i)$ with transformed weights $w_i^r = w_i^{GL} \cdot 2R/(1-x_i)^2$
2. Angular: Lebedev quadrature on the unit sphere (50–770 points, integrating spherical harmonics up to $l=11$–$29$)
3. Combined: $n_r \times n_\Omega$ points per atom: $\mathbf{r}_{ij} = r_i \hat{\Omega}_j + \mathbf{R}A$, $w{ij} = w_i^r \cdot w_j^\Omega \cdot r_i^2$
4. Molecular partitioning (Becke weights): For each atom pair, compute confocal elliptical coordinate $\mu_{AB} = (|\mathbf{r}-\mathbf{R}_A| - |\mathbf{r}-\mathbf{R}_B|)/|\mathbf{R}_A - \mathbf{R}_B|$, then smooth step function $s(\mu) = \tfrac{1}{2}[1-f(f(f(\mu)))]$ where $f(x) = (3x-x^3)/2$. The triple application of $f$ ensures smooth partition of unity.

Test 2.1 — Uniform Grid Properties ($L=10$ Bohr, $N=21$)

Test 2.2 — Uniform Grid Integration ($81^3 = 531{,}441$ pts)

Machine precision.

Test 2.3 — Becke Grid Properties ($n_\text{rad}=50$, $n_\text{ang}=194$)

The 20-order-of-magnitude weight range reflects the $r^2 dr$ volume element: tiny near the nucleus, huge far away (but integrand is zero there).

Test 2.4 — Becke Integration Accuracy (target: $\int (\alpha/\pi)^{3/2} e^{-\alpha r^2} d^3\mathbf{r} = 1$)

$10^{-14}$ accuracy regardless of exponent.

Test 2.5 — H₂ Molecule Grid

Test 2.6 — Grid Efficiency Comparison (target: error $< 10^{-5}$)

Uniform grid:

Becke grid:

Becke ~3,880 pts vs Uniform ~35,937 pts → 9.3× fewer points. The anomalous 0.30 Bohr result (worse than 0.50) is an aliasing artifact — the spacing systematically misses the Gaussian's features.

Evaluates basis function products $\chi_\mu(\mathbf{r})\chi_\nu(\mathbf{r})$ on numerical grids for overlap integrals and electron densities.

Single product: Evaluate both GTOs on all grid points, return element-wise product.

Batch shell-pair products — vectorized via NumPy broadcasting:

1. Evaluate shell $i$ → shape $(n_i, N_\text{pts})$
2. Evaluate shell $j$ → shape $(n_j, N_\text{pts})$
3. products = χi[:, newaxis, :] × χj[newaxis, :, :] → shape $(n_i, n_j, N_\text{pts})$

All $n_i \times n_j$ products computed simultaneously with no Python loops.

Overlap integral: $S_{AB} = \sum_k w_k \cdot \chi_A(\mathbf{r}_k) \cdot \chi_B(\mathbf{r}_k)$ — dot product of weighted function values. Analytical reference via the Gaussian product theorem:

$$S_{AB} = N_A N_B \left(\frac{\pi}{\alpha+\beta}\right)^{3/2} \exp!\left(-\frac{\alpha\beta}{\alpha+\beta}|\mathbf{A}-\mathbf{B}|^2\right)$$

Test 3.1 — Self-Overlap

Test 3.2 — H₂ Overlap Integral ($\mathbf{R}_A=(0,0,-0.7)$, $\mathbf{R}_B=(0,0,+0.7)$ Bohr, $\alpha=1.0$)

$S=0.375$ means 37.5% orbital overlap — characteristic of a covalent bond.

Test 3.3 — s-p Orthogonality

Zero by symmetry — $x \cdot e^{-\alpha r^2}$ is odd.

Test 3.4 — p-Shell Overlap Matrix

Max deviation from identity: $3.63\times10^{-14}$.

Test 3.5 — STO-3G Hydrogen Self-Overlap

Test 3.6 — Uniform vs Becke for Density Integration

Uniform:

Becke:

Implements three Fourier transform methods: standard FFT on uniform grids, NUFFT Type-3 via FINUFFT on Becke grids, and direct summation for validation.

FFT: Reshape density to $(N,N,N)$ → 3D FFT → scale by $h^3$ → k-grid from fftfreq(N, h) × 2π.

NUFFT Type-3 — computes $\tilde{f}(\mathbf{k}_j) = \sum_i c_i e^{i\mathbf{k}_j \cdot \mathbf{r}_i}$ where ${\mathbf{r}_i}$ are Becke points and ${\mathbf{k}_j}$ are spherical k-grid targets:

1. Weight coefficients: $c_i = f(\mathbf{r}_i) \times w_i$
2. Threshold filter: keep only $|c_i| > 10^{-14} \times \max|c|$
3. Cast to complex128 coefficients, float64 coordinates (C-contiguous via ascontiguousarray)
4. Call finufft.nufft3d3(x, y, z, c, kx, ky, kz, eps=ε, isign=1)
5. Fallback to direct summation ($O(MN)$ loop) if FINUFFT fails

FINUFFT achieves $O(N\log N + M)$ via: Gaussian gridding → FFT on oversampled grid → interpolation.

Spherical k-grid — Gauss-Legendre in each spherical coordinate, points: $\mathbf{k}{ijk} = k_i(\sin\phi_k\cos\theta_j, \sin\phi_k\sin\theta_j, \cos\phi_k)$, weights: $w{ijk} = w_i^r w_j^\theta w_k^\phi k_i^2$.

Analytical validation transforms:

$$\tilde{\chi}_s(\mathbf{k}) = N\left(\frac{\pi}{\alpha}\right)^{3/2} e^{-k^2/(4\alpha)} \cdot e^{i\mathbf{k}\cdot\mathbf{R}}, \qquad \tilde{\rho}(\mathbf{k}) = N^2\left(\frac{\pi}{2\alpha}\right)^{3/2} e^{-k^2/(8\alpha)}$$

Test 4.1 — Spherical k-Grid ($10\times10\times10$, $k_\text{max}=10$)

Test 4.2 — NUFFT vs Direct Summation

Max difference: $3.49\times10^{-13}$.

Test 4.3 — NUFFT on Normalized Gaussian ($\tilde{f}(\mathbf{k}) = e^{-k^2/(4\alpha)}$)

Test 4.4 — NUFFT on s-type GTO

Test 4.5 — Off-Center GTO Phase (center at $(1.0, 0.5, -0.3)$ Bohr)

Phase factor $e^{i\mathbf{k}\cdot\mathbf{R}}$ correctly captured.

Test 4.6 — Fourier Transform of Density ($\rho = |\chi_s|^2$, effective exponent $2\alpha$)

Test 4.7 — Convergence

$10^{-12}$ accuracy at ~35,000 points.

Computes ERIs by combining Fourier-transformed charge densities with the $4\pi/k^2$ Coulomb kernel. Implements FFT, NUFFT, and Direct methods, validated against the analytical Boys function:

$$(ab|cd) = N_a N_b N_c N_d \cdot \frac{2\pi^{5/2}}{\zeta\eta\sqrt{\zeta+\eta}} \cdot K_{ab} K_{cd} \cdot F_0(T), \qquad F_0(T) = \frac{\sqrt{\pi}}{2\sqrt{T}}\text{erf}(\sqrt{T})$$

ERI via NUFFT:

1. Weight densities: $c = \rho \odot w_\text{grid}$
2. NUFFT Type-3 transform both bra and ket densities to spherical k-grid
3. Coulomb kernel with regularization: $k^2_\text{safe} = \max(k^2, 10^{-20})$
4. Integrate: $\text{ERI} = \text{Re}!\left(\sum w_k \cdot \tilde{\rho}\text{bra} \cdot \tilde{\rho}\text{ket}^* / k^2_\text{safe}\right) / (2\pi^2)$

ERI via FFT:

1. Reshape densities to $(N,N,N)$, compute 3D FFTs scaled by $h^3$
2. k=0 singularity: set $k^2[0,0,0]=1$, compute $4\pi/k^2$, then set $\text{Coulomb}[0,0,0]=0$ (charge neutrality)
3. Sum integrand and scale by $\Delta k^3/(2\pi)^3$

ERI via Direct: Double loop over grid pairs with $r_{12}^\text{safe} = \max(|\mathbf{r}_1 - \mathbf{r}2|, 10^{-10})$ — fails catastrophically because $1/r\text{min} = 10^{10}$ corrupts the integral.

Test 5.1 — Same-Center (aa|aa)

11 significant digits.

Test 5.2 — H₂ Molecule ERIs ($\mathbf{R}_A=(0,0,-0.7)$, $\mathbf{R}_B=(0,0,+0.7)$; Grid: 69,440 pts)

Test 5.3 — Fourier vs Direct

Fourier: 16.6× faster, ~$10^{14}$× more accurate. Direct: FAILED.

Test 5.4 — Different Exponents

All ≥ 7 significant digits.

Test 5.5 — Convergence (real-space grid, fixed k-grid 3,375 pts)

Convergence (k-space grid, fixed real-space 34,720 pts)

Systematic head-to-head benchmarking of all three methods (FFT, NUFFT, Direct) on identical test cases against the analytical Boys function reference.

Test 6.1 — Same-Center (ss|ss) (Analytical: $1.1283791671$ Ha)

FFT (uniform grid):

NUFFT (Becke grid):

Direct:

CATASTROPHIC FAILURE — Coulomb singularity produces error $> 10^7$.

Test 6.2 — H₂ Molecule (bond length 1.4 Bohr)

Test 6.3 — Direct Comparison

Test 6.4 — Convergence

FFT with box size ($N=96$ fixed):

FFT error saturates at ~8% due to periodic boundary artifacts and cubic k-grid inefficiency for the spherical $1/k^2$ kernel.

NUFFT with grid refinement:

NUFFT succeeds — Fourier methods analytically handle the $1/r_{12}$ singularity through the well-behaved $4\pi/k^2$ kernel, while Becke grids concentrate every point where electron density matters.

FFT fails to converge — (1) periodic boundary artifacts, (2) finite box truncation, (3) cubic k-grid poorly samples the spherically symmetric $1/k^2$ kernel.

Direct integration fails catastrophically — when two grid points coincide, $1/r_\text{min} = 1/10^{-10} = 10^{10}$, which corrupts the entire integral.

1. Szabo A, Ostlund NS. Modern Quantum Chemistry. Courier Corporation; 2012.
2. Helgaker T, Jorgensen P, Olsen J. Molecular Electronic-Structure Theory. Wiley; 2013.
3. Becke AD. A multicenter numerical integration scheme for polyatomic molecules. J. Chem. Phys. 1988;88(4):2547–53.
4. Barnett AH, Magland J, af Klinteberg L. A parallel nonuniform fast Fourier transform library based on an "exponential of semicircle" kernel. SIAM J. Sci. Comput. 2019;41(5):C479–504.

Mehdi Hassanbeigi
Email: hasanbeigimahdi25@gmail.com

© 2026 Mehdi. All Rights Reserved.

Restrictions:

• ❌ No copying, modification, or distribution of this work is permitted