GPU Accelerated Minimal Auxiliary Basis Approach TDDFT for Large Organic Molecules

Abstract: We introduce a GPU-accelerated implementation of time-dependent density functional theory with the minimal auxiliary basis approach (TDDFT-risp) in GPU4PySCF, together with large system demonstrations carried out using the Tamm--Dancoff approximation (TDA-risp). The method combines GPU-accelerated three-center integral evaluation, tensor contractions, exchange-space truncation, omission of hydrogen atoms from the auxiliary basis, and a host memory assisted Davidson solver. On the EXTEST42 benchmark set, a conservative 40 eV exchange cutoff yields excitation-energy errors relative to standard TDA of about 0.03--0.05 eV for low-lying states. For systems of 300 to 3000 atoms, we demonstrate that TDA-risp calculations of 15 low-lying excited states with $ω$B97XD/def2-SVP complete on a single A100 GPU with wall times ranging from minutes to hours. These results position GPU-TDDFT-risp as a practical route toward excited-state calculations for large organic and biomolecular systems with thousands of atoms.

• The paper demonstrates a GPU-accelerated TDDFT approach using a minimal auxiliary basis to enable excited-state calculations for systems up to 3000 atoms with errors below 0.05 eV.
• It employs on-the-fly Coulomb evaluation, exchange-space truncation, and hydrogen exclusion strategies to significantly reduce computational cost.
• The method achieves dramatic speedups over CPU-based approaches, facilitating large-scale electronic structure and spectroscopic analyses of organic and biomolecular systems.

GPU-Accelerated Minimal Auxiliary Basis TDDFT for Large Organic Molecules

Introduction and Context

The exponential growth in computational demand for excited-state electronic structure calculations in biomolecules and organic materials, especially at scales of hundreds to thousands of atoms, has prompted the development of methods that bridge the gap between accuracy and scalability. This work presents a highly optimized implementation of linear-response time-dependent density functional theory (TDDFT) using the minimal auxiliary basis approach (TDDFT-risp), targeting commodity GPU hardware for large organic molecules (2603.29257). The method is implemented within GPU4PySCF and leverages the Tamm–Dancoff approximation (TDA-risp) for further computational efficiency.

The approach introduces algorithmic and practical advances—on-the-fly integral evaluation, memory-efficient tensor contractions, aggressive exchange-space truncation, and device/host hybrid memory strategies—to make hybrid-functional excited-state calculations feasible for systems encompassing thousands of atoms, which are otherwise infeasible with traditional CPU-based wavefunction or DFT methods.

Formalism and Algorithmic Innovations

TDDFT-risp leverages a minimal auxiliary basis composed of $s$- and $p$-type primitive Gaussian functions. This yields a reduced scaling in the evaluation of three- and four-center electron repulsion integrals (ERIs), particularly when RI techniques are employed for Coulomb and exchange terms. The TDA-risp formulation (neglecting the coupling term $\mathbf{Y}$) results in a Hermitian eigenproblem amenable to iterative solution via memory-aware Davidson algorithms.

Key algorithmic enhancements enabling large-scale practical calculations include:

• On-the-Fly Coulomb Evaluation: J-terms are constructed directly in the atomic orbital (AO) basis instead of the molecular orbital (MO) basis, avoiding storage of $\mathcal{O}(N^3)$ MO-basis tensors and exploiting AO sparsity for effective $\mathcal{O}(N^2)$–$\mathcal{O}(N^3)$ scaling.
• Exchange-Space Truncation: Exchange construction, which dominates cost for hybrid and range-separated hybrid (RSH) functionals, is truncated using orbital energy windows, significantly reducing the effective prefactor while retaining low errors.
• Hydrogen Exclusion from the Auxiliary Basis: Since hydrogens typically contribute minimally to low-lying excitations, auxiliary functions on hydrogen centers are excluded from the exchange fit, yielding substantial additional speedup for large organics.
• Memory-Aware Davidson Iterative Solver: Host/device streaming strategies enable the solution of eigenproblems whose tensors vastly exceed single-GPU memory. Davidson subspace vectors are dynamically managed in single precision, with restarts judiciously triggered before device or host memory exhaustion.

The structure of ERI tensors and exchange/Coulomb intermediates is visualized in Figure 1.

Figure 1: Schematic structure of MO-basis three-center ERI tensors for Coulomb (orange) and exchange (green) partitioning, illustrating the effect of auxiliary and MO-space truncation.

The overall computational flow, mapped to host and device memory regions and highlighting intensive data movement and contractions, is depicted in Figure 2.

Figure 2: Davidson diagonalization workflow, with major memory-resident quantities (yellow) and GPU-intensive computations/data transfers (red).

Numerical Accuracy and Validation

Benchmark Error Analysis

Rigorous accuracy benchmarks are performed on the EXTEST42 set, spanning medium-to-large organic molecules. The primary modifiable parameter is the exchange truncation window. For PBE0 and $\omega$B97XD/def2-TZVP, even a moderately aggressive 40 eV window yields excitation energy RMSE below 0.05 eV relative to full TDA Figure 3.

Figure 3: RMSE of 20 lowest excitation energies as a function of exchange window threshold for PBE0 and $\omega$B97XD.

Further, RMSE of the 20 lowest excitation energies remains below 0.05 eV for molecules of size $\geq$40 atoms—even without hydrogen auxiliary functions Figure 4.

Figure 4: RMSE of the 20 lowest excitation energies (TDA-risp vs. standard TDA) versus molecular size for the EXTEST42 set.

Systematic errors across functionals and exchange cutoffs are summarized in Figure 5.

Figure 5: Mean absolute error and RMSE for S$_1$ excitation energies, contrasting default and aggressive truncation settings, with/without hydrogen auxiliary functions.

These results validate the exchange truncation and minimal auxiliary basis for production accuracy, even in the aggressive setting (16 eV window, no H), especially for larger systems.

Spectral and Electronic Analysis

Simulated UV-vis Figure 6 and ECD Figure 7 spectra confirm that both conservative and aggressive approximation stacks preserve not just averaged state energies but also overall spectroscopic features for large molecules, with deviations for smaller systems diminishing with increasing system size.

Figure 6: UV-vis absorption spectra for the EXTEST42 set; both conservative and aggressive (16 eV/no H) TDA-risp parameterizations closely match standard TDA references.

Figure 7: Corresponding ECD spectra for EXTEST42 affirm TDA-risp's fidelity for rotational strength calculations.

Charge-transfer (CT) excitations, frequently sensitive to long-range exchange, are also faithfully captured: in prototypical Betaine 30 and D-$p$0-A systems, the spatial hole-electron separation and dominant NTOs are preserved under the risp approximations Figure 8.

Figure 8: NTO and hole–electron density analyses for charge-transfer states in Betaine 30 and D-$p$1-A systems demonstrate accurate long-range coupling.

Performance and Scalability

Timing and Resource Analysis

Wall-time profiling shows that 15-state TDA-risp calculations on $p$2–$p$3 atom systems are routinely achievable on a single NVIDIA A100 in wall times from minutes to a few hours—orders of magnitude faster than conventional CPU-based RIJCOSX approaches.

The relative speedup compared to both state-of-the-art CPU (PySCF) and established quantum chemistry packages (ORCA) is dramatic, as in Figure 9.

Figure 9: Wall-time speedup for TDA-risp (GPU and CPU) over RIJCOSX TDA (ORCA, 32 core MPI), peaking at $p$4 for 480-atom systems.

Even for systems at the upper end of EXTEST42 (480 atoms), TDA-risp closely tracks standard RIJCOSX spectra and S$p$5 energies within 0.02 eV Figure 10.

Figure 10: UV absorption spectrum of a 480-atom tetramer, showing near-identity between TDA-risp and conventional ORCA RIJCOSX spectra.

Large-System Davidson Solver Behavior

Davidson iteration histories for large aggregates and proteins exhibit complex convergence patterns, often requiring tens of iterations, multiple restarts, and careful residual monitoring for robust eigenvalue solution Figure 11.

Figure 11: Residual norms of the Davidson solver as a function of iteration for large aggregates and proteins, illustrating convergence complexity and the effect of restarts.

Large-System Excited-State Analysis

The ability to perform full-system TDDFT excited-state analysis on chromophores embedded in realistic environments is demonstrated on systems such as fluorescent protein crystal structures (PDB IDs 5EXB, 5EXC, $p$63145 atoms) and the Photosystem II reaction center (1299 atoms). Hole–electron analyses display local and CT-like excitations spanning the protein and chromophore domains without any fragmentation approximation or implicit environmental modeling Figure 12.

Figure 12: Hole–electron density decompositions for S$p$7 of fluorescent protein 5EXC and S$p$8 of a Photosystem II model, resolved at the thousand-atom scale.

Conclusion

This work establishes that with a minimal auxiliary basis, AO-domain contraction, strategic exchange truncation, and GPU-centric implementation, TDDFT and TDA can be rendered practical for excited-state calculations in large organic and biomolecular systems at scales (up to $p$93000 atoms) inaccessible to previous first-principles methodologies. Errors relative to standard TDA are below 0.05 eV under conservative settings, and speedups versus CPU implementations exceed two orders of magnitude for large systems.

The main computational bottleneck for even larger systems ($\mathbf{Y}$03000 atoms) shifts to the SCF ground-state step due to hardware and software constraints in ground-state diagonalizers. The theoretical and practical impact is wide-ranging: screening of red-shifted/broadened chromophores, simulation of organic device aggregates, and entire-protein quantum photophysics. Immediate future directions include robust preconditioning, extension to open-shell and radical systems, analytical gradients, automated auxiliary basis adaptation, and multi-GPU distribution for systems exceeding ten thousand atoms.