FourPhonon_GPU: Accelerated Phonon Transport

• The paper demonstrates FourPhonon_GPU achieving >25x acceleration in four-phonon scattering rate evaluation, with thermal conductivity accuracy within 0.1%.
• It employs a heterogeneous CPU–GPU strategy using OpenACC kernels to overcome the N^4 scaling challenge in dense Brillouin-zone meshes.
• The framework retains the established FourPhonon physics under the RTA, enabling accurate, high-fidelity calculations for complex material systems.

FourPhonon_GPU is a GPU-accelerated extension of the FourPhonon package for computing three-phonon (3ph) and four-phonon (4ph) scattering rates and lattice thermal conductivity from first principles. It targets the principal computational bottleneck in phonon transport modeling: the combinatorial growth of scattering processes as the Brillouin-zone mesh and the number of phonon branches increase. For a $16\times 16\times 16$ silicon mesh, the reported process counts are on the order of $9.0\times 10^5$ for three-phonon scattering and $7.6\times 10^9$ for four-phonon scattering, a regime in which exhaustive CPU evaluation can consume thousands of core-hours. FourPhonon_GPU preserves the validated physics and interfaces of FourPhonon while using a heterogeneous CPU–GPU strategy and OpenACC kernels to accelerate the numerically intensive stages without loss of accuracy, reporting more than $25\times$ acceleration in scattering-rate kernels, more than $10\times$ end-to-end runtime reduction, and agreement with the CPU baseline within $0.1\%$ in thermal conductivity (Guo et al., 1 Oct 2025).

1. Computational role and problem setting

FourPhonon_GPU addresses the prohibitive cost of four-phonon scattering, whose process count scales as $N^4$ with the number of q-points $N$. In conventional CPU/MPI workflows, the main difficulties are energy and momentum filtering across very large candidate spaces, nested loops over branches and q-points, and matrix-element evaluation for a vast number of allowed processes. The framework is designed to preserve the full physics of the underlying FourPhonon formalism: it does not introduce sampling or machine-learning surrogates, and it explicitly accelerates mode-resolved 3ph and 4ph scattering rate evaluation under the relaxation-time approximation (RTA), weighted phase space (WP) integrals, and end-to-end conductivity calculations including both 3ph and 4ph contributions (Guo et al., 1 Oct 2025).

Its significance lies in the specific status of four-phonon scattering in thermal transport. The reported motivation is that including 4ph scattering is essential for accurate $\kappa$ in many materials at room and higher temperatures, but the computational burden had made dense-mesh and complex-crystal calculations difficult in practice. The framework is therefore positioned not as a new physical approximation, but as a performance-oriented implementation that retains FourPhonon’s established formulations while reorganizing the computation to match GPU hardware.

2. Scattering formalism and numerical treatment

The framework evaluates phonon scattering rates on a discrete q-mesh under energy and momentum conservation. For three-phonon processes, the reported schematic expression is

$$\Gamma^{(3)}_{\mathbf{q}s} = \frac{2\pi}{\hbar} \sum_{\mathbf{q}'s',\mathbf{q}''s''} \left| V^{(3)}(\mathbf{q}s,\mathbf{q}'s',\mathbf{q}''s'') \right|^2 \Big[ (n_{\mathbf{q}'s'}+n_{\mathbf{q}''s''}+1)\, \delta(\omega_{\mathbf{q}s}-\omega_{\mathbf{q}'s'}-\omega_{\mathbf{q}''s''}) + (n_{\mathbf{q}''s''}-n_{\mathbf{q}'s'})\, \delta(\omega_{\mathbf{q}s}+\omega_{\mathbf{q}'s'}-\omega_{\mathbf{q}''s''}) \Big],$$

with momentum conservation $9.0\times 10^5$0. For four-phonon processes, the reported general structure is

$9.0\times 10^5$1

with momentum conservation $9.0\times 10^5$2 (Guo et al., 1 Oct 2025).

Delta functions are treated numerically by Gaussian broadening. The width $9.0\times 10^5$3 is computed from differences in group velocities, and the paper states that a broadening factor of $9.0\times 10^5$4 is used when 4ph is included. The weighted phase space combines Bose factors, matrix elements, and broadened delta functions, and is computed alongside the scattering rates. Thermal conductivity is then assembled under the RTA from mode heat capacities, group velocities, and lifetimes, with lifetimes determined by the sum of three-phonon, four-phonon, and optional additional scattering rates such as isotopes and boundaries.

The matrix elements retain the standard lattice-dynamical dependence on interatomic force constants (IFCs). Third-order IFCs determine $9.0\times 10^5$5, while fourth-order IFCs determine $9.0\times 10^5$6. The implementation is stated to follow the FourPhonon formulations of Feng and Ruan for four-phonon channels and selection rules. In this sense, FourPhonon_GPU changes the execution model rather than the underlying phonon-scattering physics.

3. Heterogeneous CPU–GPU execution model

The central algorithmic choice is to remove branching-heavy and irregular operations from the GPU. The reported core idea is to eliminate GPU warp divergence by moving all filtering with substantial branching to the CPU, and to send to the GPU only valid, energy- and symmetry-allowed events. CPU tasks therefore include enumeration of symmetry-allowed candidate tuples satisfying momentum conservation, energy checks using Gaussian broadening, and construction of compact index arrays such as Ind2_all, Ind3_all, analogous 4ph structures, cumulative offsets N_accum, and per-mode counts N_plus. Harmonic properties, phase-space handling, and post-processing that remain inherently CPU-based in FourPhonon are also kept on the host (Guo et al., 1 Oct 2025).

GPU tasks are restricted to massively parallel work: evaluation of weighted phase-space contributions, Bose factors, and matrix elements for allowed processes, followed by reductions into mode-resolved $9.0\times 10^5$7 and WP arrays. The reported rationale is that the decoupled nature of scattering processes under RTA matches the SIMT model, so millions of independent threads can compute per-process contributions efficiently once invalid processes have already been removed. A recurring practical point in the benchmarks is that GPU-only filtering is not generally advantageous; the hybrid CPU–GPU path is specifically reported to avoid the severe warp divergence that appears for sparse 4ph processes.

OpenACC is the implementation vehicle. The paper describes directive-based acceleration in Fortran loops to minimize code restructuring and preserve MPI compatibility. Nested loops over $9.0\times 10^5$8, branch indices, and the additional index required for 4ph are flattened using collapse clauses; loop order is rearranged to promote coalesced access to frequencies, velocities, and eigenvectors; and computations for $9.0\times 10^5$9 and matrix elements are inlined to avoid function-call overhead. Reductions for $7.6\times 10^9$0 and WP use OpenACC reductions rather than atomics.

4. Data structures, scheduling modes, and memory behavior

The data layout is built around compact representations of allowed processes. Ind2_all and Ind3_all store flattened indices, while N_accum provides prefix sums for mode-local segments. This organization allows contiguous storage of accepted processes per mode and therefore more coalesced thread traversal on the device (Guo et al., 1 Oct 2025).

Two GPU scheduling modes are provided. In the all-modes strategy, all data and indices are preloaded and a single kernel is launched across all modes and processes. This offers the highest throughput, but also the highest memory footprint. In the mode-by-mode strategy, per-mode indices and data are copied to the GPU and a kernel is launched for each mode. This reduces memory demand but increases host–device transfers and launch overhead. The paper reports that all-modes can exceed available memory for dense meshes; for example, $7.6\times 10^9$1ph+$7.6\times 10^9$2ph at $7.6\times 10^9$3 exceeded $7.6\times 10^9$4 GB, and all-modes on a $7.6\times 10^9$5 mesh can require approximately $7.6\times 10^9$6 GB.

The implementation is synchronous: CPU enumeration completes before GPU kernels begin. Preloading of constant data such as frequencies, velocities, eigenvectors, and IFC-derived data is used to reduce PCIe traffic. In all-modes, the indices are copied once and reused; in mode-wise execution, only the per-mode segments are transferred. The paper identifies asynchrony and overlap as future work rather than a current capability.

A further implementation detail concerns symmetry-restricted iteration spaces. Four-phonon combinations can occupy triangular iteration regions, whereas OpenACC loop collapsing requires rectangular domains. The adopted solution is to perform triangular-set enumeration on the CPU and pass only valid indices to the GPU. This is a specific example of the broader design principle: regularize the work before it reaches the accelerator.

5. Benchmark results and hardware dependence

The reported benchmarks focus on silicon at $7.6\times 10^9$7 K under the RTA and show that the acceleration does not come at the expense of scientific agreement. Thermal conductivity deviates by less than $7.6\times 10^9$8 from the CPU-only FourPhonon baseline, with minor variation attributed to compiler differences (Guo et al., 1 Oct 2025).

Quantity	Reported result
End-to-end runtime	$7.6\times 10^9$9 faster for both 3ph and 3ph+4ph across tested meshes
Scattering kernels	$25\times$0 speedup for 3ph and $25\times$1 speedup for 4ph
All-modes vs mode-by-mode	Additional $25\times$2 speedup for 3ph at $25\times$3; $25\times$4 for 3ph+4ph at $25\times$5
Hybrid CPU–GPU vs GPU-only	$25\times$6 faster for 3ph at $25\times$7; GPU-only degraded below CPU baseline for 3ph+4ph at $25\times$8
Accuracy	$25\times$9 deviation in $10\times$0 compared to CPU-only FourPhonon

The architecture dependence is explicitly linked to double-precision capability. The paper reports kernel performance ordering as A100 $10\times$1 $10\times$2 A30 $10\times$3 $10\times$4 A10 $10\times$5, and states that FP64 capability dominates. This is consistent with the use of double precision throughout. A common misconception in GPU computing is that a larger number of GPU threads is sufficient to guarantee speedup; the reported results instead emphasize that filtering strategy, divergence control, memory capacity, and FP64 throughput are decisive for this workload.

6. Inputs, outputs, limitations, and research significance

FourPhonon_GPU is implemented inside FourPhonon and preserves original inputs, outputs, and file formats. Inputs include harmonic phonon dispersions and group velocities from second-order IFCs, obtained through DFPT or finite displacement, third- and fourth-order IFCs for matrix elements, symmetry operations, q-point mesh specification, temperatures, and optional scattering sources such as isotopes and boundaries when total rates require them. FourPhonon routines perform harmonic-property and symmetry-grid preprocessing, while FourPhonon_GPU adds CPU enumeration of allowed 3ph and 4ph processes under energy and momentum constraints. Outputs include mode-resolved $10\times$6, $10\times$7, total scattering rates, weighted phase-space arrays, and the thermal conductivity tensor versus temperature under the RTA (Guo et al., 1 Oct 2025).

The practical workflow is reported as follows: harmonic properties are precomputed on the CPU; allowed processes are enumerated on the CPU and recorded in contiguous arrays; constant arrays are preloaded to the GPU; OpenACC kernels compute Bose factors, WP contributions, and inlined matrix elements from IFC-derived data; and the CPU then performs unit conversion, mesh averaging, summation of 3ph and 4ph rates, lifetime construction, and assembly of $10\times$8. For deployment, the paper recommends the heterogeneous option as the default, mode-by-mode execution as the initial choice for large problems, and all-modes execution when memory permits. GPU runs typically use one CPU core plus one GPU. The reported software environment includes the NVIDIA HPC compiler with OpenACC support for GPU builds, Intel/OneAPI or GNU for CPU-only builds, and benchmark builds using CUDA 12.6 and nvc 23.5.

The principal limitations are also explicit. The present implementation is restricted to the RTA; iterative BTE solvers, including variational or relaxation-time improved schemes, are described as future work because of memory and host–device transfer constraints. All-modes execution can exceed $10\times$9–$0.1\%$0 GB for dense meshes with 4ph, so mode-wise fallback is necessary in some cases. FP64 is treated as required for accuracy, though mixed precision is identified as a possible future direction. Broadening magnitude affects both physics and performance; the paper states that $0.1\%$1 would yield even higher kernel acceleration, but must be physically justified. Finally, the quality of fourth-order IFCs remains sensitive to DFT choices and exchange-correlation functionals, which is a general FourPhonon consideration rather than a GPU-specific issue.

In broader computational-physics terms, FourPhonon_GPU is important because it makes fully resolved four-phonon transport calculations routine on modern HPC systems without introducing approximations beyond those already present in the RTA workflow. The reported implications are larger meshes, more branches, faster temperature- and material-space studies, and high-fidelity baselines for validating surrogates and assessing higher-order IFC sensitivity. A plausible implication is that its main contribution is methodological alignment between phonon-scattering algorithms and accelerator hardware: control-heavy selection remains on the CPU, while embarrassingly parallel numerical kernels are moved to the GPU.