CFFT Acceleration

Overview

Standard FFT-based SCFT solvers operate on the full $N^3$ grid, which contains redundant information when the crystal phase has high space group symmetry. Polyorder.jl integrates CrystallographicFFT.jl to exploit this redundancy through a Crystallographic FFT (CFFT) mode, operating on a symmetry-reduced subgrid $M^3$ where $M_i = N_i / L_i$.

Key benefits:

1. Universal Symmetry Support: Works for all 230 space groups — including those with glide/screw symmetry (e.g., Gyroid SG 230 $Ia\bar{3}d$) that DCT cannot handle.
2. Significant Speedup: MDE solving runs ~$L^3 \times$ faster. For Pm-3m ($L=2$): ~8×; for Ia-3d ($L=2$): ~8× with additional star compression. A further ~2× from real FFT (rfft/irfft) vs complex FFT.
3. Memory Reduction: Solver buffers reduced by $L^3 \times$. Combined with SymmetricStorage: up to 99% total memory savings.
4. Machine Precision: Mathematically equivalent to full FFT solver ($|\Delta F| < 10^{-14}$).

Mathematical Basis

Subgrid Decomposition

For a space group $G$ with point group $P$ acting on an $N^3$ grid, the symmetry-independent information lives on a subgrid of size $M^3 = (N/L)^3$, where $L$ is the stride factor determined by the space group and grid dimensions.

The CFFT algorithm (KRFFT) decomposes the full $N^3$ DFT into:

1. An $M^3$ FFT on the subgrid
2. A reconstruction step that maps $M^3$ spectral coefficients to the $n_\text{spec}$ points of the spectral ASU (Asymmetric Unit)

This is exact: the subgrid retains all symmetry-independent spatial information, and reconstruction recovers the full spectrum at uniquely determined ASU wave vectors.

Spectral ASU Operations

The spectral ASU is the minimal set of independent Fourier modes under the space group. All spectral operations (diffusion kernel $e^{-k^2 \Delta s}$, stress tensor) can be computed directly on the spectral ASU without ever touching the full $N^3$ grid:

• Diffusion kernel (cfft_k2): $k^2$ values at spectral ASU points
• Stress tensor (cfft_kk_orbsum): Orbit-averaged $k_\alpha k_\beta$ for each Voigt component

MDE Solving on Subgrid

The Operator Splitting Formula (OSF) for the Modified Diffusion Equation:

\[q(\mathbf{r}, s + \Delta s) = e^{-w \Delta s / 2} \cdot \text{ICFFT}\left[ e^{-k^2 b^2 \Delta s / 6} \cdot \text{CFFT}\left[ e^{-w \Delta s / 2} \cdot q(\mathbf{r}, s) \right] \right]\]

All operations occur on the $M^3$ subgrid (real-space multiplication) and $n_\text{spec}$-point spectral ASU (spectral multiplication). The full $N^3$ grid is never allocated or accessed.

Architecture

CFFTAuxiliaryField

The core type representing a field on the subgrid:

struct CFFTAuxiliaryField{...} <: AbstractAuxiliaryField
    data::A                  # Size: M³ = (Nx/L, Ny/L, Nz/L). CPU Array or CuArray.
    lattice::BravaisLattice  # Full lattice (with physical dimensions)
    cfft_pair::PP            # Bidirectional real CFFT plan (rcfft! / ircfft!)
    sg_num::Int              # Space group number (1–230)
    N_full::NTuple{3, Int}   # Full grid dimensions
    star_map::SubgridStarMap # Weighted integration map
end

All internal buffers (stress tensor, spectral work arrays, propagator storage) are device-aware: they are allocated on the same device as data using KernelAbstractions.allocate or similar(data, ...), ensuring seamless CPU/GPU operation.

Pipeline: MDE Solve

q(M³) ──→ exp(-w·ds/2) (pointwise, M³)
       ──→ rcfft! (M³ → n_spec complex, via rfft)
       ──→ exp(-k²·ds) (pointwise, n_spec)
       ──→ ircfft! (n_spec → M³, via irfft)
       ──→ exp(-w·ds/2) (pointwise, M³)
──→ q(M³, s+ds)

Pipeline: Stress Tensor (Spectral ASU)

Instead of FFT on $N^3$, stress is computed entirely in the spectral ASU:

\[S_{\alpha\beta}(s) = -\frac{1}{N^6} \sum_{i=1}^{n_\text{spec}} \overline{\hat{q}_i(s)} \cdot \hat{q}^\dagger_i(s) \cdot m_i \cdot \overline{k_\alpha k_\beta}_i\]

where $m_i$ is the orbit multiplicity and $\overline{k_\alpha k_\beta}_i$ is the orbit-averaged $k_\alpha k_\beta$ from cfft_kk_orbsum. The inner reduction uses GPU-compatible broadcasting (sum(real.(conj.(F̂_q) .* F̂_qc) .* mult .* kk_j)) to avoid scalar indexing. CFFTStressData stores kk_orbsum and multiplicity as device-adapted vectors via Adapt.adapt.

Weighted Statistics

Observables computed on the subgrid require a volume weight correction because the subgrid density of points is $L^3$ times lower than the full grid:

mean(w::CFFTAuxiliaryField) = mean(w.data) * (N_vol / M_vol)^(1/D)
# N_vol = prod(N_full), M_vol = prod(size(data))

This ensures Hw, Hs, F, and residuals are exactly consistent with full-grid FFT values.

Spectral Updater Support (SIS/ETD/PO)

All SpectralSCFTAlgorithm updaters now support CFFT fields natively via trait-based dispatch:

The mechanism uses a spectral_transform_type trait that returns Val(:fft) for standard fields and Val(:cfft) for CFFTAuxiliaryField. Each updater's step! dispatches to the appropriate implementation:

• FFT path (default): Uses plan_fft/plan_ifft on $N^3$ grid
• CFFT path: Uses rcfft!/ircfft! on $n_\text{spec}$-point spectral ASU

The κ response coefficients (from RPA) are computed on the spectral ASU using cfft_k2 — this is mathematically equivalent since κ depends only on $k^2$ and the spectral ASU orbits represent unique $k^2$ values.

Nested updaters (Anderson, NGMRES) that use spectral updaters as preconditioners work transparently — they operate in real space on the subgrid and delegate spectral operations to the preconditioner's CFFT path.

Integration with Performance Profiles

CFFT is compatible with all performance profiles that use symmetry:

CFFT and DCT are mutually exclusive — the constructor automatically selects CFFT when a CFFTAuxiliaryField is provided, bypassing DCT detection.

Verification

• Free Energy: CFFT matches FFT to machine precision ($|\Delta F| < 10^{-14}$)
• Stress Tensor: CFFT stress matches FFT stress to $4.3 \times 10^{-19}$ (machine limit)
• Density: Pointwise difference $< 10^{-15}$
• GPU: Validated on NVIDIA RTX 2080 Ti ($|\Delta F_\text{GPU-CPU}| < 10^{-11}$)
• GPU VariableCell: Full GPU support including stress tensor, resize, and grid hysteresis
• Spectral Updaters: CFFT SIS/ETD/ETDPEC/PO/SISF/ETDF produce identical convergence to FFT counterparts

Supported Space Groups

CFFT supports all 230 space groups in 3D. The speedup depends on the stride factor $L$:

Performance Benchmarks

Benchmarked on an AB diblock copolymer system with SG 230 (Ia-3d Gyroid), ds=0.01, profile=:balanced. The benchmark measures update_propagator! — one full MDE propagation step.

CPU (single-threaded)

GPU (NVIDIA RTX 2080 Ti)

Initialization Overhead

CFFT initialization (planning, star map construction) is faster than standard FFT initialization after CrystallographicFFT.jl v0.2.1 optimization:

Benchmark scripts: benchmark/cfft_performance.jl (CPU), benchmark/cfft_performance_gpu.jl (GPU), benchmark/cfft_init_profile.jl (init profiling).

Comparison with DCT Acceleration