EOM-CC Excited-State Gradients and Nonadiabatic Couplings on a Consumer GPU from a Contraction-DAG with Laplace-Transform J/K Kernels

Abstract: We present a unified, memory-bounded GPU realization of equation-of-motion coupled-cluster (EOM-CC) excited-state gradients and interstate nonadiabatic couplings (NACMEs) on a single 8\,GB consumer GPU. Both are built from one contraction directed acyclic graph: the EOM-CC relaxation is the reverse-mode transpose of the forward density build rather than a per-state re-derivation, and an atomic-orbital-direct Laplace-transform $J/K$ kernel, made non-symmetric ($J^x(A,B)\neq J<sup>x(B,A)$) by the transition densities, resolves every energy denominator with no four-index molecular-orbital tensor; a two-sided Davidson returns both eigenvectors from one device-resident, spin-pure solve. The pipeline is \emph{validated end to end at small scale}: gradients and NACMEs match finite differences across four spin multiplicities and full configuration interaction to $&lt;!10<sup>{-12}$ for two electrons, and the excited-state gradient matches the independent \textsc{Psi4} code to $\le!4.6\times10<sup>{-7}~E_h/a_0$ from \ce{H2O} to aromatic benzene. The kernels and the ground-state solve reach chromophores ($\le!730$ AO) in 8\,GB, and a frozen-natural-virtual compression lets the eigensolver \emph{execute} a complete excited-state gradient and $Q$--$B$ NACME of the chlorophyll-core chromophore \ce{Mg}-porphine (def2-SVP, $439$ AO) on the card. We present that run as a \emph{capability demonstration} -- executed and translationally invariant to machine zero, but anchored only piece-wise and bounded by a direct convergence study at ${\sim}10<sup>{-2}~E_h/a_0$ -- not a converged spectroscopic result. The validated small-scale capability and the memory-bounded implementation are the contribution.

• The paper presents a GPU-optimized framework that integrates analytic EOM-CCSD gradients and NACMEs via a contraction-DAG and Laplace-transform based J/K kernel.
• It employs reverse-mode automatic differentiation and Cholesky decomposition to achieve high accuracy and efficiency within an 8 GB memory envelope.
• Numerical results validate speedups up to 43× over CPU baselines and confirm method consistency with established quantum chemical techniques.

Unified EOM-CCSD Gradients and Nonadiabatic Couplings on a Consumer GPU: Contraction-DAG with Laplace-Transform J/K Kernels

Overview

This work presents an integrated, memory-bounded GPU implementation of equation-of-motion coupled cluster (EOM-CCSD) excited-state analytic gradients and interstate nonadiabatic coupling matrix elements (NACMEs), delivering both properties from a single contraction-directed acyclic graph (DAG) with a Laplace-transform-based J/K kernel. All computations, including demanding correlated excited-state response properties, fit within the 8 GB memory envelope of a consumer-level NVIDIA GPU. The architecture synthesizes established theoretical constructs—automatic differentiation of contraction graphs, AO-direct Laplace quadrature for energy denominators, and non-symmetric transition-density contractions—into an end-to-end, device-resident solution accessible on commodity hardware.

Theoretical Framework

EOM-CC Excited-State Gradients and NACMEs

EOM-CC methods compute excitation energies and response properties by representing excited states with linear excitation operators acting on a correlated reference. The key challenge addressed is evaluating analytic nuclear gradients and interstate NACMEs, which require both state-specific and transition density matrices, and careful treatment of non-Hermitian left/right eigenvector structures.

The gradient functional contracts relaxed one- and two-particle densities with derivative integrals and J/K builds; every term's orbital-energy denominator is resolved via a Laplace transform, reformulating conventional four-index molecular-orbital (MO) tensor operations into sequences of atomic-orbital (AO) contractions over modest quadrature grids.

Laplace-Transform AO-Direct J/K Kernels

A central component is the AO-direct Laplace J/K kernel. Integration over energy denominators is rewritten as finite exponential sums using minimax quadrature, and the resulting AO-density contractions are handled efficiently on the GPU, avoiding the memory bottleneck of storing $O(N^4)$ MO tensors. This approach systematically reduces the memory scaling at each correlation level.

Contraction DAG and Relaxation as Transpose

All tensor contractions involved in building densities and their response are compiled into a contraction DAG. The analytic gradient and NACME machinery exploit reverse-mode automatic differentiation: the relaxation equations (amplitude-response, Z-vector) are generated as the DAG's reverse-mode transpose. This approach unifies the computation of gradients and NACMEs within one abstraction, exploiting efficient kernel reuse and automatic path optimization.

Figure 2: The HGP-OS recurrence kernel: each CUDA thread evaluates a full ERI quartet DAG, fully asynchronously and linearly scalable, with regime-specific handling for angular momenta.

Figure 4: Coulomb (J) and exchange (K) energy builds execute as cuBLAS contraction DAGs on a Cholesky factor $B$, with the N$^4$ ERI tensor never formed; $J$ and $K$ operations remain bandwidth-bound.

Non-symmetric Transition-Density Builds

Key to general EOM-CC response theory is the necessity of handling non-symmetric transition densities ($J^x(A,B) \neq J^x(B,A)$), stemming from biorthonormal left/right states. This requirement is met by designing ordered J/K kernels accepting two independent density arguments, both for the forward (density) and backward (relaxation) passes.

Device Implementation Strategy

Kernel and Memory Optimization

All kernel operations are tailored to the memory and bandwidth constraints of a consumer-grade GPU. The contraction-DAG-driven chunking partitions kernel execution to fit tight register and shared-memory budgets, minimizing spills and ensuring high SM occupancy. The AO integrals are evaluated via an in-house HGP-OS code generator rather than preexisting (e.g. TeraChem) engines, allowing regime-specific unrolling, tiling, or looped execution tailored by angular momentum.

Figure 6: Gradient execution: forward and energy-weighted density assembly DAGs run concurrently on the GPU, feeding derivative-ERI kernels for analytic gradients.

Cholesky Decomposition and CuBLAS Integration

The dominant two-electron integral contractions are handled efficiently via Cholesky decomposition of the ERI tensor, reducing them to matrix multiplications that execute on cuBLAS. Only the necessary three-index Cholesky factors $B$ and two-index densities are retained in memory, never the full four-index ERI tensor, respecting the strict device footprint limit.

Unified CPU/GPU Validation Stack

The GPU path is validated in three complementary ways: agreement with dense CPU implementations on small molecules, analytic vs. finite-difference gradient/NACME agreement across multiple spin multiplicities, and exact correspondence with FCI for two-electron systems. At production scale, kernel equivalence across implementations is used for correctness, along with stringent checks on translation invariance and internal consistency.

Numerical Results and Performance

Accuracy Validation

Maximum analytic–finite-difference deviations are $\lesssim 10^{-6}~E_h/a_0$ for gradients and NACMEs across molecules spanning all spin multiplicities, including singlet, doublet, triplet, and quartet cases and cross-validated against external codes (Psi4). Analytic gradients coincide with FCI up to $<10^{-12}$ $E_h/a_0$ for two-electron systems.

Chromophore-Scale Capability Demonstration

The pipeline is exercised in production for full analytic excited-state gradients and $B$0–$B$1 NACMEs of Mg-porphine (def2-SVP, 439 AO), entirely within 8 GB. The final calculation includes complete per-atom gradient vectors and NACME components, with force completeness certified by machine-precision translation invariance. Because external validation is infeasible at this size, method-intrinsic uncertainties due to virtual space truncation are explicitly bounded.

Figure 8: Kernel roofline for the RTX 4060: J/K and gradient kernels achieve high bandwidth efficiency but J remains bandwidth-limited and is retained in FP64.

Figure 10: Cost decomposition for chromophore Born–Oppenheimer MD steps shows $B$2-density contractions dominate analytic gradient runtime, and two-electron terms dominate both energy and gradient builds.

Figure 1: Per-SCF-iteration breakdown by energy term and walltime across benchmark systems up to 730 AO.

Figure 3: Scaling analysis reveals near-quadratic J and cubic K walltime dependence on AO count, with gradient assembly maintaining favorable scaling under memory and bandwidth constraints.

Figure 5: Analytic two-electron gradient assembly achieves $B$3–$B$4 speedups over single-thread CPU baseline for feasible systems.

Figure 7: In-method timing for transition-density response gradient kernel matches ground-state kernel to within 1%, confirming kernel-level unification and efficiency.

Kernel Throughput and Bottlenecks

Achieved bandwidth utilization for J/K and $B$5-density build kernels approaches 80–94% of the theoretical DRAM roof, with the main bottleneck being the inherently bandwidth-limited nature of the AO-direct J build. The three-center gradient term and Cholesky metric eigenproblem are the remaining cost drivers.

Memory Management and Trajectory Readiness

Geometry-frozen Cholesky metric and screening caches yield $B$6 per-step speedup for excited-state MD applications, reusing device-resident intermediates across multiple gradients and NACME evaluations at fixed geometry. All device-resident intermediates are ephemeral and released immediately after their utility, maintaining a sawtooth memory pattern peaking strictly below the card limit.

Implications and Future Perspectives

The presented device pipeline makes high-level correlated excited-state nuclear gradients and nonadiabatic couplings, previously restricted to datacenter hardware, available on affordable consumer GPUs. The contraction-DAG spine, Laplace AO kernels, and DAG-transpose relaxation are designed to generalize, by construction, across MPn, CCSDTQ, and (structurally, though not demonstrated) multireference (CASSCF, CASPT2) families.

Practical implications include enabling full quantum dynamics and photochemistry simulations for nontrivial chromophores on previously unavailable hardware; methodologically, the unification of kernel signature for symmetric and non-symmetric/transition density builds streamlines code generation, maintenance, and extensibility. Future work includes extending the approach to cavity-QED polaritonic Hamiltonians and possibly adapting relaxation-as-transpose to RPA or Green's functions regimes, contingent on underlying self-consistency structures.

Conclusion

The work establishes and validates a robust, memory-bounded GPU pipeline for EOM-CCSD excited-state gradients and interstate nonadiabatic couplings. Utilizing a single contraction graph abstraction, reverse-mode automatic differentiation, and AO-direct Laplace J/K kernels, the approach attains high accuracy, stringent validations, and substantial speedups for realistic systems on commodity hardware. The implications are both practical—democratization of quantum photochemical simulation capabilities—and conceptual, in enabling modular, graph-theoretic automatic differentiation strategies for correlated quantum response properties.