EOM-CC Excited-State Gradients and Nonadiabatic Couplings on a Consumer GPU from a Contraction-DAG with Laplace-Transform arXiv:2607.05622v1 [physics.chem-ph] 6 Jul 2026
J/K Kernels Ruben Dario Guerrero∗,†,‡ †NeuroTechNet S.A.S., 1108831, Bogotá, Colombia ‡Quantum and Computational Chemistry Group (QCCG), Universidad Nacional de Colombia, Bogotá, Colombia E-mail: [email protected]
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) ̸= J x (B, A)) by the transition densities, resolves every energy denominator with no four-index molecularorbital tensor; a two-sided Davidson returns both eigenvectors from one device-resident, spin-pure solve. The pipeline is validated end to end at small scale: gradients and NACMEs match finite differences across four spin multiplicities and full configuration interaction to < 10−12 for two electrons, and the excited-state gradient matches the
1
independent Psi4 code to ≤ 4.6 × 10−7 Eh /a0 from H2 O to aromatic benzene. The kernels and the ground-state solve reach chromophores (≤ 730 AO) in 8 GB, and a frozennatural-virtual compression lets the eigensolver execute a complete excited-state gradient and Q–B NACME of the chlorophyll-core chromophore Mg-porphine (def2-SVP, 439 AO) on the card. We present that run as a capability demonstration—executed and translationally invariant to machine zero, but anchored only piece-wise and bounded by a direct convergence study at ∼10−2 Eh /a0 —not a converged spectroscopic result. The validated small-scale capability and the memory-bounded implementation are the contribution.
1
Introduction
When a molecule absorbs light it evolves on excited-state potential-energy surfaces and can return to the ground state through conical intersections — geometries where two electronic states become degenerate and the Born–Oppenheimer separation breaks down. Simulating this photochemistry — vision, photosynthesis, photovoltaics, DNA photoprotection — requires two quantities at every nuclear geometry: the gradient of each excited state (the force that drives the dynamics) and the nonadiabatic coupling matrix element (NACME, the rate at which neighbouring states exchange population near an intersection). Computing both accurately for correlated wavefunctions is the central bottleneck of nonadiabatic dynamics. With equation-of-motion coupled cluster (EOM-CC), 1,2 the benchmark single-reference method for excited states, obtaining these quantities has traditionally meant re-deriving the response (relaxation) equations by hand for each state and running them on datacenter GPUs. Here we obtain both automatically — as the reverse-mode transpose (the same mechanism as backpropagation in machine learning) of a single contraction graph — and engineer the resulting non-symmetric Coulomb/exchange kernels to run within the 8 GB budget of a consumer NVIDIA GeForce RTX 4060 (GPU performance is reported in Sect. 4.4), putting a capability that today needs a cluster allocation onto a card a student already owns. 2
Each ingredient of this construction has a mature literature, and we engage it directly so that the contribution is not overstated. EOM-CC excited-state gradients and nonadiabatic couplings are established quantities: the analytic EOM-CC gradient and its Lagrangian (Z-vector) formulation were developed by Stanton and Gauss 3,4 and implemented for spinconserving and spin-flip states by Krylov and co-workers, 5 while EOM-CC nonadiabatic and derivative couplings have been formulated and validated against multireference references by the Krylov and Koch groups. 6–9 The reverse-mode view of relaxation is likewise not new: differentiable programming through electronic structure now yields nuclear gradients, response densities, and derivative couplings by automatic differentiation, from the Hartree–Fock proof of concept of Tamayo-Mendoza et al. 10 to the PySCFAD framework of Zhang and Chan; 11 in particular, Zhang et al. 12 obtain coupled-cluster response properties directly from the reverse pass in place of an explicit Z-vector solve — the ground-state realization of “relaxation = reverse pass.” Compiling many-body equations into a condensed, order-optimized contraction graph traces to the Tensor Contraction Engine 13 and continues in modern coupled-cluster code generators, and resolving correlated denominators directly in the atomic-orbital basis by a Laplace transform underlies the AO-direct Laplace MP2 gradient of Ochsenfeld and co-workers. 14 What we add is the specialization of these threads to the non-Hermitian EOM-CC response problem: a single graph transpose that emits both the right-state amplitude response (excited-state gradients) and the interstate left/right Zvector (NACMEs), uniformly across EE/IP/EA rather than per state and per property, fed through an AO-direct minimax-Laplace J/K build made non-symmetric by the transition densities (J x (A, B) ̸= J x (B, A)), and delivered within the 8 GB of a consumer GPU. To our knowledge this combination is without direct precedent; the individual mechanisms are not. The gap is therefore not in any one ingredient but in their union for the non-Hermitian excited-state-response problem. Correlated nonadiabatic dynamics through conical intersections needs excited-state gradients and interstate NACMEs together, yet no existing route to them combines all three of: (i) the relaxation obtained as a single reverse-mode
3
graph transpose, uniform across EE/IP/EA and across gradients and NACMEs, rather than per-state, per-property hand-derived Lagrangian/Z-vector machinery; (ii) a fully AOdirect correlated pipeline, with the non-symmetric transition-density Coulomb/exchange and their energy denominators resolved in the atomic-orbital basis, never forming a four-index molecular-orbital tensor; and (iii) a memory footprint that fits the 8 GB of a consumer GPU rather than workstation or datacenter hardware. The adjacent threads each leave this union open: differentiable-programming electronic structure 10–12 supplies reverse-mode response for ground-state and mean-field methods on general frameworks, not the non-Hermitian EOM-CC response and not AO-direct under a memory budget; the established EOM-CC gradient and NACME implementations 3–9 work in the molecular-orbital basis, per state and per property, and target capable hardware; AO-direct Laplace correlation 14 addresses ground-state MP2 gradients, not excited-state transition densities; and contraction-graph code generation 13 emits equations and kernels, not the response-as-transpose nor the ordered, non-symmetric transition-density J/K. The consequence is that correlated excitedstate forces and couplings — the inputs nonadiabatic dynamics consumes — remain practically out of reach on the hardware most groups own. Closing that gap is the point of this work. Our thesis is that the correlated excited-state forces and couplings above can be brought within the memory budget of a single consumer GPU, and that one contraction DAG is the mechanism that makes it so: its reverse-mode transpose is the EOM-CC relaxation, so a single abstraction yields both the excited-state gradient and NACME (as the DAG transpose) and the spill-bounded CUDA J/K kernels for symmetric (ground-state) and non-symmetric (transition/interstate) densities alike. The contribution is threefold, ordered by what is new. (1) The systems capability: a device-native, AO-direct, non-symmetric J/K kernel (J x (A, B) ̸= J x (B, A), never forming a four-index molecular-orbital tensor) that realizes the complete non-Hermitian EOM-CCSD excited-state gradient and interstate NACME within the 8 GB of a consumer GPU, its DAG-
4
derived chunking reducing register/shared-memory spilling (e.g. the (dd|dd) live set from 96 to 34 KB) so that the memory-bound transition-density builds run compute/bandwidthbound (Sects. 2.4, 3). (2) The enabling mechanism: a single contraction DAG whose reversemode transpose is the EOM-CC relaxation—specializing the coupled-cluster differentiation identity to the non-Hermitian EOM case—from which the symmetric, non-symmetric, and gradient builds are all emitted (Sect. 2). (3) The validation stack: finite-difference and determinant-oracle checks across all four spin multiplicities, an FD-free exact full-CI crosscheck, a ζ=0 control, and an independent cross-code (Psi4) check reaching the excited-state gradient itself from H2 O to aromatic benzene (Sect. 4.3), with commodity-GPU roofline and throughput at chromophore scale (Sects. 3–4). We are explicit about what is validated versus what is demonstrated. The validated contribution is twofold: the memory-bounded kernel architecture, and the end-to-end correctness of the integrated gradient and NACME at small scale—where finite differences across four multiplicities, an exact full-CI identity, a determinant oracle, and an independent cross-code check (against Psi4, reaching the excitedstate gradient itself and holding from H2 O to aromatic benzene) all apply. The complete chromophore-scale excited-state gradient and Q–B interstate NACME of Mg-porphine (Table 4) are then executed within 8 GB—the novel systems capability, a calculation that today needs a cluster allocation—but, no external oracle being affordable at that size, are presented as a bounded capability demonstration (with the ∼10−2 Eh /a0 FNO-truncation uncertainty of Sect. 4.3 stated throughout), not a spectroscopically converged result. The LT-AO J/K kernel and the ground-state CD-RCCSD solve independently reach 730 AO on the card. Because canonical EOM-CC implementations cannot run EOM-CCSD at chromophore scale on commodity hardware, no external large-scale comparison exists; at that scale the method is supported by finite differences, internal consistency, and the kernel-identity argument together with the small-scale and benzene-scale cross-code agreement above. Section 2 develops the theory—the EOM-CC gradient functional (Sect. 2.1), the LT-AO J/K kernel that resolves every energy denominator (Sect. 2.2), the relaxation as the DAG transpose
5
(Sect. 2.3), and the non-symmetric build the transition densities force (Sect. 2.4); Sect. 3 gives the computational details, and Sect. 4 the validation and GPU results.
2
Theory
The theory is organized around a single observation: every analytic first derivative in this framework is the contraction of one- and two-particle (transition) density matrices with skeleton (derivative) integrals through one Coulomb/exchange (J/K) kernel, with every energy denominator resolved by a Laplace quadrature into a short grid of atomic-orbital (AO) density builds. We develop this in four steps. Section 2.1 defines the EOM-CC excitedstate densities and the gradient functional they enter. Section 2.2 develops the AO-direct Laplace-transform (LT-AO) J/K kernel that resolves every denominator and assembles the gradient—the formalism that keeps the calculation within an 8 GB memory budget. Section 2.3 shows that the relaxation (response) equations are the transpose of the forward density-build graph, generated mechanically rather than re-derived per state. Section 2.4 treats the non-symmetric J/K build that the EOM transition densities force.
Notation and conventions. Indices i, j, k, l label spin orbitals occupied in the Hartree– Fock reference |0⟩ ≡ |Φ0 ⟩, a, b, c, d the unoccupied (virtual) orbitals, p, q, r, s either, and µ, ν, κ, ρ the atomic-orbital (AO) basis; Cµp are the molecular-orbital (MO) coefficients and εp the orbital energies. We work in second quantization, with p̂† and p̂ the creation and annihilation operators for spin orbital p, and hpq and (pq|rs) (chemist notation) the oneand two-electron MO integrals. Excitations are generated by the operators
τ̂ia = ↠î,
τ̂ijab = ↠b̂† ĵ î,
...
(1)
(each creates particles in virtual orbitals and holes in occupied ones); the full set is written {τ̂µ }, and |µ⟩ ≡ τ̂µ |0⟩ are the excited determinants. The cluster operator is T̂ = 6
1 a a ia ti τ̂i + 4
P
P
ab ab ijab tij τ̂ij + · · · , with amplitudes t ≡ {tµ } fixed by the coupled-cluster equa-
tions ⟨µ|H̄|0⟩ = 0; the EOM operators R̂k , L̂k (Eq. (2)) are expanded in the same manifold {τ̂µ }, with r0 the weight on |0⟩. A superscript x denotes a skeleton derivative with respect to a nuclear coordinate x at fixed MO coefficients — the core-Hamiltonian derivative hxµν , x the overlap derivative Sµν , and the integral derivative (µν|κρ)x — while orbital relaxation
is restored separately through the coupled-perturbed Hartree–Fock (CPHF) coefficients U x (Eq. (8)). One- and two-particle density matrices are γpq = ⟨p̂† q̂⟩ and Γpqrs = ⟨p̂† r̂† ŝ q̂⟩; their interstate (transition) forms are defined in Eq. (6). Convention for the density contractions and J/K builds. The wavefunction operators above are written in spin orbitals, but the one- and two-particle density contractions, the cumulant collapse, and the Coulomb/exchange (J/K) builds and gradient assembly that the GPU kernels implement (Eqs. (7), (8), (19), (20), (23), and (24)) are stated in the spatial closed-shell (RCCSD) convention that the production code uses—in which the separable twoparticle density carries the spatial exchange coefficient 12 and the closed-shell two-electron build is 2J − K on the per-spatial-orbital density, equivalently J − 12 K on the spin-summed total density (the two related by 2J − K = 2(J − 12 K)). Throughout, J x and K x denote the single bare builds of Eq. (18), which carry no embedded prefactor; the exchange coefficient 1 , the closed-shell spin sum (the half-vs-total-density factor of two), and the bra–ket pair2
interchange (pq|rs) = (rs|pq) are written explicitly wherever they occur and are never folded into the build. The spin-orbital expressions enter only the reference verification engine (Sect. 4.1), which uses a separate spin-orbital implementation.
2.1
EOM-CC excited states, response densities, and the gradient functional
We begin with the object being differentiated: the EOM-CC excited-state energy, its left and right eigenvectors, and the non-symmetric densities that the gradient contracts.
7
Setup. EOM-CC parameterizes an excited state by acting on the coupled-cluster ground P ab ab P τ̂ij + · · · , 1 state with a linear, CI-like excitation operator R̂k = r0 + ia ria τ̂ia + 41 rij |Ψ0 ⟩ = eT̂ |Φ0 ⟩,
|Ψk ⟩ = R̂k |Ψ0 ⟩,
(2)
where |Ψ0 ⟩ is the CC ground state built from the reference |Φ0 ⟩. Inserting this ansatz into the Schrödinger equation gives Ĥ eT̂ R̂k |Φ0 ⟩ = Ek eT̂ R̂k |Φ0 ⟩, with Ek the total energy of state k. Because R̂k and T̂ are both pure excitation operators they commute, [T̂ , R̂k ] = 0, so R̂k slides through eT̂ and the latter factors out. Left-multiplying by e−T̂ then yields an eigenvalue problem, e−T̂ Ĥ eT̂ − Ek R̂k |0⟩ = 0,
(3)
in which the bare Hamiltonian is replaced by the similarity-transformed operator H̄ ≡ e−T̂ Ĥ eT̂ , whose eigenvalue Ek = E0 + ωk is the total energy of state k (E0 the coupledcluster ground-state energy, the R̂0 = 1 root; ωk the excitation energy). This is the central object of EOM-CC: a similarity transform leaves the spectrum unchanged, so the construction is formally exact; in practice T̂ and R̂k are truncated to a common excitation level (singles through quadruples here, EOM-CCSDTQ). The price of the transform is that H̄ is not Hermitian—the wave operator eT̂ is nonunitary. Its left and right eigenvectors are therefore distinct objects (unlike ordinary quantum mechanics, where one bra/ket pair suffices), and they form a biorthonormal pair ⟨Lj |Rk ⟩ = δjk rather than an orthonormal one. That non-Hermiticity, and the L ̸= R structure it forces, is the root of nearly everything “non-symmetric” in this paper. Collecting the twosided eigenproblem,
H̄ = e−T̂ Ĥ eT̂ ,
H̄ R̂k |0⟩ = (E0 + ωk )R̂k |0⟩,
⟨0|L̂k H̄ = (E0 + ωk )⟨0|L̂k ,
(4)
biorthonormalized as ⟨Lj |Rk ⟩ = δjk . In the k-fold excitation basis {τ̂µ }, H̄ is represented by
8
the (non-symmetric) EOM Jacobian
Aµν = ⟨µ| [ H̄, τ̂ν ] |0⟩,
(5)
whose right and left eigenvectors are Rk , Lk with eigenvalue ωk . A is the operator that the relaxation step below transposes. For EOM-CCSDTQ the operators R̂k , L̂k span the singles through quadruples excitation manifolds (Fig. 1).
R̂1 : singles
R̂3 : triples
R̂2 : doubles
R̂4 : quadruples
Figure 1: EOM-CCSDTQ excitation manifolds. left) eigenvectors of H̄ span P a aThe1 right P ab(and ab singles through quadruples, R̂k = r0 + ia ri τ̂i + 4 rij τ̂ij + · · · ; the dashed line is the Fermi level, filled circles are reference electrons, and red arrows promote n electrons from occupied to virtual orbitals. The non-symmetric transition densities γ AB , ΓAB of Eq. (6) are formed from these manifolds. The quadruples manifold is the reach of the symbolic reference-equation generator (p† q) and the determinant oracle; the device pipeline and every validated gradient and NACME in this work are EOM-CCSD.
Why the densities are non-symmetric. Because H̄ (hence A) is non-Hermitian, the left and right eigenvectors differ (Lk ̸= Rk ), so the one- and two-particle (transition) densities AB γpq = ⟨LA | p† q |RB ⟩,
† † ΓAB pqrs = ⟨LA | p r s q |RB ⟩
(6)
are non-symmetric: γ AB ̸= (γ AB )⊤ . Operationally this forbids the usual A ↔ B permutation symmetry of the Coulomb/exchange build, forcing the ordered kernel J x (A, B) ̸= J x (B, A) of Sect. 2.4 below — the GPU kernel at the heart of this work. The common shortcut Lk = Rk symmetrizes the densities and corrupts the gradient; we instead retain the true biorthonormal triple (ωk , Rk , Lk ). We obtain it from a native two-sided non-Hermitian Davidson that advances the right and left Krylov subspaces of H̄ together and returns both eigenvectors 9
from a single solve. Because the pair is non-Hermitian, ordinary Gram–Schmidt fails: we use biorthogonal Gram–Schmidt (deflating each Lk against the converged right set), holding the right vector at unit L2 norm and rescaling only the dual left vector so that ⟨Lk |Rk ⟩ = 1; this fixes the relative normalization of the biorthonormal pair on which every transition density (and hence the gradient) depends. The solver is device-resident: each matvec A r and its transpose A⊤ w are evaluated by the same GPU kernels that assemble the gradient, so one solve yields (ωk , Rk , Lk ) directly on the card. This eliminates the two-step seed—a right vector from a borrowed dense Davidson followed by a GMRES inverse-iteration left vector—used to bootstrap earlier validation: on H2 O/STO-3G the native pair reproduces that borrowed pair to ∆ω = 4.5 × 10−9 , with right- and left-vector overlaps of unity and biorthonormality error ∥L⊤ R − I∥ at 10−16 . The solver is multiplicity-general. For a closed-shell reference it targets a chosen total spin Q through a matrix-free spin-sector projector PS = S ′ ̸=S (Ŝ 2 − σ ′ )/(σ − σ ′ ) with σ = S(S + 1); because [H̄, Ŝ 2 ] = 0 the projection commutes through the iteration and costs no additional H̄ matvec, separating singlets from triplets in a single run. For an open-shell reference the spin-orbital singles-and-doubles manifold is not closed under Ŝ 2 , 15 so no determinant-space projector can reach a spin-pure state; we instead iterate in a spin-adapted (genealogical) configuration-state-function (CSF) basis on a semicanonical ROHF reference, which renders doublets and quartets spin-pure by construction (Sect. 4.2). The validated device path is exercised at small scale—H2 O/STO-3G for the closed-shell device gates, and gate-scale (nso ≤ 18) open-shell systems for the multiplicity-general path (Table 3). At chromophore scale the full-virtual (uncompressed) eigensolve is a demonstrated capability rather than a tolerance-converged production result; the FNO-compressed production eigensolve used here converges to a < 5 meV drift—looser than the 2.3 × 10−6 Eh benzene/STO-3G class, but two orders of magnitude below the 1.5 eV Q–B gap and well inside the separation of each bright 1 Eu state from its neighbours, so the state assignment—the observable used here—is robust even where the absolute excitation energy is not tolerance-converged (Sect. 4.4).
10
Gradient as a density functional. The excitation-energy gradient is the derivative of the density functional
F AB [f, g; t] =
X
AB fpq γpq + 12
pq
X pqrs
gpqrs ΓAB pqrs − sAB Ecorr ,
(7)
where fpq = hpq and gpqrs = (pq|rs) are the one- and two-electron Hamiltonian integrals, ΓAB is the two-particle transition density of Eq. (6), Ecorr = ⟨0|H̄|0⟩ − EHF is the groundstate correlation energy, and sAB = ⟨LA |RB ⟩ is the interstate overlap (sAB = 1 for a state gradient, 0 for a symmetry-distinct NACME). The last term removes the disconnected reference correlation that the cluster operator injects into the densities; in the normal-ordered representation the subtracted pieces are t1 in the ov block of γ AB and the connected doubles combination −(tijab + tia tjb − tib tja ) in the oovv block of ΓAB . At the converged cluster amplitudes t this makes F AB = ωk in value (the disconnected piece vanishes; verified to 3.5 × 10−9 ) while its amplitude derivative ∂F AB /∂t ̸= 0—the non-vanishing derivative is exactly what the amplitude-response equation of Sect. 2.3 below is constructed to cancel.
Master gradient and NACME. The relaxed densities feed one AO-direct contraction, X dEA X x x Wµν Sµν + J x (γ, γ) − 21 K x (γ, γ) + 2 Tr[GF U x ] + Enuc γµν hxµν + = , dx µν µν
(8)
where γ, Γ here denote the relaxed one- and two-particle densities (γ ≡ P of Eq. (9) below, and its two-particle partner), W is the energy-weighted density (conjugate to the overlap derivative S x ), GF the generalized Fock and U x the coupled-perturbed Hartree–Fock (CPHF) orbital-response coefficients (the 2 Tr[GF U x ] term is the orbital Z-vector contribution), and each energy denominator inside γ, Γ, W is resolved by a Laplace grid of quadrature nodes {τα , wα } (the τα are scalar Laplace nodes, unrelated to the excitation operators τ̂µ ; the P α wα is suppressed here for compactness). Section 2.2 expands this denominator resolution explicitly and writes Eq. (8) in its Laplace-grid form, Eq. (20). 11
The symbols entering Eqs. (8) and (20) are the relaxed densities. The total one-particle density decomposes additively,
P = DGS (T, Λ) + DEE (T, L, R) + Dζ (T, ζ),
(9)
into the ground-state (Λ-relaxed) density, the EOM transition density built from L, R, and the amplitude-response density Dζ obtained by evaluating the ordinary Λ-linear coupledcluster response density with the multiplier ζ of Sect. 2.3 below; P is thus the γ of Eq. (8) P once relaxation is included. The energy-weighted density is Wµν = − pq εp Ppq Cµp Cνq (in x ), and the Laplace form εp → εp eεp τα ; only its symmetric part contributes to the symmetric Sµν P P generalized (non-symmetric) Fock matrix is GFpq = r hpr Prq + rst (pr|st) Γqrst (free index
p on the one- and two-electron integrals, contracted against the relaxed densities), whose occupied–virtual block is the right-hand side of the orbital Z-vector (CPHF) equation that produces the orbital-response coefficients U x . Because L ̸= R† , the two-particle transition density ΓAB additionally carries the blocks γ (2) |ovoo and γ (2) |vvov that are absent from any symmetric L = R treatment. Analytic EOM-CC nonadiabatic and derivative couplings have themselves been formulated and validated against multireference references by several groups, 7,8 including the biorthonormal CCSD construction of the interstate left/right elements 9 closest to the nonHermitian build used here. The quasidiabatic interstate NACME follows by the off-diagonal Hellmann–Feynman route 6 with the transition densities γ AB , ΓAB in place of state densities,
λAB = x
⟨LA | ∂x H̄ |RB ⟩ , EB − E A
(10)
the 1/(EB − EA ) pole being the conical-intersection (Berry-connection) singularity. Equation (10) is the quasidiabatic coupling: L, R are frozen at the reference geometry, which removes the eigenvector-response terms (EA − EB ) ∂x ⟨LA |RB ⟩. For different-irrep pairs sAB = ⟨LA |RB ⟩ = 0, so the Ecorr subtraction and the nuclear-repulsion term drop out; 12
same-irrep pairs (sAB ̸= 0) carry an additional gap term that the present implementation in general; the symmetrized ̸= λBA guards by assertion rather than evaluates. Note λAB x x p BA value we report is the geometric mean λAB x λx , formed in post-processing. The device kernel returns only the one-sided coupling numerator ⟨LA |∂x H̄|RB ⟩; the 1/(EB − EA ) gap division and the geometric-mean symmetrization are applied afterward by the caller in an FP64 island, not on-device. Both gradient and NACME reduce to the same AO-direct J/K kernel of Sect. 2.2 fed different densities.
2.2
The AO-direct Laplace-transform J/K kernel
Equations (8) and (10) reduce every gradient and NACME to density contractions in which each energy denominator must be resolved. We now develop the kernel that performs that resolution entirely in the atomic-orbital (AO) basis, never forming a four-index molecularorbital (MO) tensor. This is the formalism that keeps the calculation within an 8 GB memory budget; we present it as a derivation, from the bottleneck to the master Laplace-grid gradient.
The denominator bottleneck. A correlated gradient contracts amplitudes and densities against two-electron integrals divided by orbital-energy denominators D (e.g. D = εa + εb − εi − εj for a double excitation, εa − εi for a single; ε the canonical orbital energies, i, j occupied and a, b virtual). Carried out conventionally, this forms and stores four-index MO tensors, of which the all-virtual (vvvv) block dominates memory and scales as the fourth power of the virtual-space dimension. On commodity hardware this block alone exceeds the available memory well before the molecules of interest are reached. The Laplace identity. The Häser–Almlöf identity 16,17 replaces the reciprocal denominator by an integral, approximated on a short minimax (best-L∞ exponential-sum) quadrature grid {τα , wα }, 18,19
1 = D
Z ∞ 0
e
−Dτ
dτ ≈
13
nτ X α=1
wα e−Dτα ,
(11)
with nτ of order ten points sufficing for sub-microhartree accuracy; the grid is generated for a target accuracy over the spanned denominator range [Dmin , Dmax ], over which the minimax error decays nearly exponentially in nτ , 18 controlled by nτ together with the ratio Dmax /Dmin , not by nτ alone. The decisive property is that the exponential of the (additive) denominator factorizes into one-index factors,
e−Dτα = e+εi τα e+εj τα e−εa τα e−εb τα ,
(12)
so each orbital index carries its own scalar weight. The contraction can then run index by index in the AO basis, and no four-index MO tensor is ever built.
Laplace-scaled coefficients and the AO density. The per-index factors of Eq. (12) are absorbed directly into the MO coefficients, defining Laplace-scaled occupied and virtual coefficients at each grid point α,
α Ṽµa = Cµa e−εa τα .
α Ũµi = Cµi e+εi τα ,
(13)
Contracting these scaled coefficients over the orbital indices yields an AO-basis density at grid point α, (1;α) Dµν =
X
α α Ũµi Ṽνa .
(14)
ia
Equation (14) is the level-1 density. For correlation order n the construction iterates to a depth L = L(n) (L = n − 1 for n = 2, 3, L = ⌊n/2⌋ for n ≥ 4), one Laplace grid per (k)
denominator, the level-k density carrying the order-k amplitude residual Ria between the scaled coefficients, (k;α1 ,...,αk ) Dµν =
X ia
(k;α ,...,αk−1 )
αk Ũµi Ria 1
αk Ṽνa ,
(15)
(k)
with Ria the occupied–virtual projection of the residual produced by the forward J/K cascade through the preceding grid points (Eqs. (17)–(19)). The ket that enters the gradient
14
adds to this forward density the backward density D̂(L) constructed as the DAG transpose (Sect. 2.3 below), Ξ(L;α) = D(L;α) + D̂(L;α) ,
(16)
which carries the relaxation (Λ / Z-vector / ζ) contribution. The forward cascade (Eqs. (13)– (15)) is lines 5–10 of Algorithm 1.
The J and K builds. Each AO density is contracted with the two-electron integrals through the Coulomb and exchange builds
J(D)µν =
X κρ
(µν|κρ) Dκρ ,
X
K(D)µν =
κρ
(µκ|νρ) Dκρ ,
(17)
which are the standard AO-direct Coulomb/exchange operations 20 evaluated with GPU Gaussian-integral engines. 21 Resolving a correlated denominator directly in the AO basis through such builds, with no AO→MO transform, follows the AO-direct Laplace MP2 gradient of Schweizer, Doser, and Ochsenfeld; 14 the present kernel extends that pseudodensity construction to the non-Hermitian transition densities below. The gradient requires the AOdirect derivative (skeleton) builds, in which the integral derivatives (µν|κρ)x are contracted against two (generally distinct) density arguments,
J x (A, B) =
X (µν|κρ)x Aµν Bκρ ,
K x (A, B) =
X (µκ|νρ)x Aµν Bκρ .
(18)
µνκρ
µνκρ
These derivative builds are the J/K kernel of this work; here A, B denote the two AO density arguments (not the EOM Jacobian A of Eq. (5)). Equation (18) is the ordered build whose non-symmetry, J x (A, B) ̸= J x (B, A), the EOM transition densities force (Sect. 2.4 below). The undifferentiated forms of Eq. (17) additionally enter the generalized (Laplace) Fock back-transform, α F̃pq =
X µν
Cµp 2J(D(1;α) ) − K(D(1;α) ) µν Cνq ,
15
(19)
which returns each grid contribution to the MO basis for the response steps of Sect. 2.3 below.
The master LT-AO gradient. Assembling the per-grid builds gives the explicit Laplacegrid form of the master gradient, Eq. (8), in its relaxed AO-direct realization, h i X dEA X Y x (1;α1 ) (L;α) x (1;α1 ) (L;α) x + Pµν hxµν +Wµν Sµν , (20) ,Ξ −K D ,Ξ = wαk 2J D dx µν α k where P and W are the relaxed one-particle and energy-weighted densities and Ξ(L;α) is the combined forward+backward AO density of Eq. (16). The per-grid build is written here as 2J x − K x on the per-spatial (half) densities D(1;α1 ) , Ξ(L;α) ; this is the same closed-shell two-electron build as the J x − 12 K x of the master gradient Eq. (8) on the spin-summed total density (2J − K = 2(J − 21 K), the factor 2 being the spin sum), with J x , K x the bare builds of Eq. (18). The multi-index α = (α1 , . . . ) runs over the Laplace grids of all denominators, Q and k wαk collects their quadrature weights. Cost. Each reciprocal denominator 1/D is thus replaced by a short sum of J/K builds, one per grid point, each of cost O(N 2 ) in the AO dimension N for a direct (integral-recomputing) scheme, the asymptotic cost once Schwarz screening is applied (the formal cost of one undifferentiated build is O(N 4 )). The four-index vvvv block is never formed: because each build (Eq. (17)) contracts the integrals against a two-index AO density, the all-virtual MO summation is never assembled—it is replaced by the AO-direct build. At the MP2 level the conventional O(N 5 ) denominator contraction becomes nτ AO-direct J/K builds; at the doubles level and above, the all-virtual block that dominates conventional memory is eliminated outright. This replacement of stored four-index tensors by recomputed two-index AO builds is what makes the excited-state gradients of this work fit within an 8 GB budget.
16
EOM-CC specialization. The excited-state denominators of EOM-CC carry the excitation energy ωk , 1/(Dk − ωk ), where Dk is the orbital-energy denominator of the excitation class (a single εa − εi , a double εa + εb − εi − εj , etc.). The same quadrature applies with a single, ω-dependent reweighting of the grid, X 1 ≈ w eα (ω) e−Dk τα , Dk − ω α
w eα (ω) = wα eωτα ,
(21)
valid for ω < mink Dk , i.e. the regime set by the HOMO–LUMO gap; for higher (Rydberg/chargetransfer) roots approaching ω → Dk a contour shift in the complex-τ plane keeps the integral finite, or the IP-EOM treatment applies. Only the per-grid weight changes: the AO-direct J/K kernel of Eqs. (17)–(18) and the grid points τα are unchanged. The EOM-CC excitedstate gradient therefore uses the same LT-AO kernel as the ground state, fed the transition densities of Eq. (6) in place of the state densities. The procedure, end to end. Algorithm 1 collects the complete gradient/NACME procedure— the forward density cascade of this section, the reverse-mode relaxation derived in Sect. 2.3 below, and their assembly through the J/K kernel of Sect. 2.4—into one numbered pass, with each step cross-referenced to its working equation. The remaining two subsections supply the content the algorithm invokes by name: why the reverse pass (lines 11–13) is guaranteed to be the transpose of the forward pass (Sect. 2.3), and why a single ordered kernel covers every build the algorithm calls (Sect. 2.4).
2.3
Relaxation as the contraction-DAG transpose
The two foregoing subsections fixed what must be contracted—the transition densities and the gradient functional (Sect. 2.1)—and how every energy denominator inside them is resolved into AO-direct J/K builds (Sect. 2.2). What remains is the relaxation: the response of the amplitudes to nuclear motion, traditionally the most laborious part of an analytic gradient. We now make the “DAG is the spine” thesis concrete and show that the relaxation 17
Algorithm 1 EOM-CC excited-state gradient via the AO-LT contraction DAG. Lines 5– 10 are the forward density cascade (Sect. 2.2); lines 11–13 are its reverse-mode transpose (Sect. 2.3); the closing builds reuse the single ordered J/K kernel (Sect. 2.4). The contraction DAG that condenses, differentiates, and orders these steps is compiled once and reused across the grid and across methods. Require: geometry; converged CC amplitudes T (⟨µ|H̄|0⟩ = 0); target root k (state A ≡ k) 1: Solve the two-sided EOM eigenproblem for the biorthonormal R̂k , L̂k ▷ Eqs. (4), (5) 2: Form the non-symmetric transition densities γ AB , ΓAB ▷ Eq. (6) 3: Subtract the disconnected Ecorr pieces (t1 ; connected doubles) ▷ Eq. (7) 4: Generate the Laplace grid {τα , wα } for the spanned denominator range ▷ Eq. (11) 5: for each grid multi-index α = (α1 , . . . , αL ) do 6: Form Laplace-scaled coefficients Ũ α , Ṽ α ▷ Eq. (13) 7: Build the level-1 AO density D(1;α1 ) ▷ Eq. (14) 8: for k = 2, . . . , L do ▷ level-k cascade (k) 9: J/K-build → Laplace–Fock back-transform → residual R ▷ Eqs. (17), (19) 10: Build the level-k AO density D(k;α) ▷ Eq. (15) 11: Reverse pass = DAG transpose: seed F̄ = 1, adjoint each forward edge → source ξ = ∂F AB /∂t ▷ Eq. (22) ⊤ (L) 12: Solve the amplitude response A ζ = −ξ; build the backward density D̂ ▷ Eq. (22) 13: Combine the ket Ξ(L;α) = D (L;α) + D̂ (L;α) ▷ Eq. (16) x 14: Assemble the relaxed P, W and orbital (U ) response ▷ Eqs. (8), (9) X x (1) dEA XY x = wαk 2J (D , Ξ)−K x (D(1) , Ξ) + Pµν hxµν +Wµν Sµν ▷ Eq. (20) 15: return dx α k µν NACME variant: run lines 1–15 with mixed (L̂A , R̂B ) densities and the 1/(EB − EA ) factor ▷ Eq. (10)
18
equations fall out as the transpose of the forward density-build graph, with no per-state re-derivation.
The five-slot template. We read any correlation method as five interchangeable parts; the first four are method-specific, the fifth universal. S1 reference and orbital partition (a single Hartree–Fock determinant for EOM-CC); S2 amplitudes and coefficients (the cluster operator T̂ and the biorthonormal left/right eigenvectors L̂k , R̂k for EOM-CC); S3 oneand two-particle densities (the non-symmetric transition densities γ AB , ΓAB for couplings); S4 relaxation/response (the self-consistent ζ amplitude-response, developed below); S5 the universal atomic-orbital Laplace-transform (LT-AO) J/K kernel of Sect. 2.2. Only S1–S4 change between methods; S5 and the contraction-DAG compiler that condenses, differentiates, and orders it are shared — which is why the same engine is designed to specialize across MP2 through CASPT2. That generality is a design property of the template; in this work we develop and validate only the EOM-CC case. The labels S1–S5 index the stages of Fig. 2. The density build and its contraction through the LT-AO kernel (Fig. 2) are never manipulated term by term: the working expressions are compiled once into a single contraction DAG GL whose nodes are intermediate tensors and whose edges are multilinear contractions. Three properties of that one object carry the leverage of this work, and each maps onto a step the algorithm above takes for granted. Condensation. The hundreds of contractions in a truncated EOM-CC density share many subexpressions — the same half-contracted integral block, the same T2 -dressed intermediate — that a literal term-by-term evaluation would rebuild repeatedly. Compiling to a DAG merges identical subexpressions into shared nodes, so each is computed once and reused by every term that needs it. The amortization compounds along the two axes the algorithm loops over: a shared node is computed once and reused across all Laplace-grid points α (the ω-reweighting of Eq. (21) touches only the per-grid scalar weight, not the node), and once a
19
method’s leaves are swapped in (S1–S4 of the template) the same condensed skeleton serves MPn, CC, and EOM builds. Transpose. The reverse of GL is the relaxation (Λ / Z-vector / ζ) equation, read off by reverse-mode without per-state re-derivation (developed below). Because the transpose is taken on the graph rather than on a hand-written term list, gradients and NACMEs — and every EE/IP/EA variant at every truncation level — inherit it as one mechanism rather than a family of derivations. Path. The contraction order is optimized once (dynamic programming over the linear extensions of GL ) and that single schedule is reused across the Laplace grid and across methods (Fig. 3), so the optimization cost is paid once and amortized over every grid point and every downstream method. Compiling a complex computation into a DAG and exploiting its graph-theoretic structure is the canonical strategy of optimizing compilers, where a basic block is represented as a DAG whose shared nodes are common subexpressions and the compiler performs common-subexpression elimination, dead-code elimination, and instruction scheduling by topological ordering of that graph. 22 Within quantum chemistry, this compile-to-graphthen-condense-and-order strategy was established for many-body methods by the Tensor Contraction Engine, 13 which compiles coupled-cluster and many-body equations into a contraction graph, factors out common subexpressions, and optimizes the contraction order — the direct ancestor of the Condensation and Path properties below. The three properties above are exactly these established operations specialized to the density build: Condensation is common-subexpression elimination, Path is instruction scheduling (a topological order of GL ), and Transpose is reverse-mode differentiation of the graph, the same graph rewrite that powers automatic-differentiation and machine-learning frameworks (the backpropagation invoked above). The same toolkit recurs across disciplines—optimizing compilers, automaticdifferentiation frameworks, and coupled-cluster code generators. 23,24 In quantum chemistry specifically, graphs already drive automated GPU code generation: the automated code
20
engine of Song, Wang, and Martı́nez represents a single integral-evaluation program as a dataflow graph and transforms it to enumerate store-versus-recompute code variants, autotuned for a target device. 25 A complementary strand applies the same idea one level down, at the recurrence relations themselves: a layered (topological-level) code generator emits optimized implementations of arbitrary recurrences—orthogonal polynomials, special functions, and the molecular-integral recurrences (McMurchie–Davidson, Obara–Saika, Rys, and the Boys function)—and is reported to exceed expert hand-optimization. 26 That recurrence-level codegen is the lineage of the per-quartet HGP-OS recurrence kernels emitted by the code generator of §3 (Fig. 6). In both cases the graph organizes the evaluation of a fixed integral kernel; it is a different object from the contraction spine we use here. Our DAG is the manybody density-contraction itself, and its graph-theoretic structure is the method—shared-node condensation, one reusable contraction path, and a reverse-mode transpose that, because the graph builds the response functional F AB , is the relaxation, generating the response equations (ζ amplitude response, interstate Z-vector) rather than re-deriving them. This machinery is inherited and mature; we claim no novelty for DAGs, common-subexpression elimination, scheduling, reverse-mode differentiation, or graph-driven GPU code generation. Our contribution is their specialization to the correlated excited-state response setting — the EOM-CC gradient and NACME — and specifically the identity that the DAG transpose is the relaxation, realized through Algorithm 1. These properties are not incidental to coupled-cluster structure. GL mirrors the structure of the Goldstone/MBPT diagrammatic expansion: each diagram is a node, each contracted internal line an edge, and a shared diagram fragment maps to a shared DAG node — so condensation is the diagrammatic factorization made automatic, the transpose is the conjugate (response) diagram set generated mechanically, and the optimized path is a contraction ordering of the diagrams. GL is the diagrammatic structure made executable, differentiable, and optimizable.
21
method-specific (S1–S4)
Contraction-DAG compiler (structural invariant): condense shared intermediates · transpose ⇒ S4 · optimize path once, reuse over grid × methods S1–S2 reference + amplitudes / coefficients: T (n) (MPn) | T, L, R (CC/EOM) | CI c / MPS (MR)
transpose ⇒ S4
S3 one- / two-particle density matrices γ, Γ (state, or transition γ AB , ΓAB ) S4 relaxation ⇒ non-symmetric γ, Γ, W (Λ/Z | ζ | Z+CP-MCSCF | none)
X α
skeleton integrals hx , S x , (µν|κρ)x
R ∞ −Dτ 1 dτ D= 0 e
Laplace
S5 UNIVERSAL AO-LT J/K KERNEL X X wα [ 2J x − K x ](D, D′ ) + hx γ + S xW gradient dEA /dx
NACME λAB x
or
Figure 2: The universal gradient/NACME pipeline. Stages S1–S4 (orange) are methodspecific—they decide which densities are built and how they are relaxed; for EOM-CC, S2 supplies the biorthonormal T, L, R, S3 the non-symmetric transition densities, and S4 the self-consistent ζ response. Stage S5 (blue) is the same kernel for every method: a Laplace-resolved, non-symmetric J/K contraction. Spanning the pipeline, the contractionDAG compiler (violet) generates the S4 relaxation equations as its transpose (dashed edge). workingequation term list (hundreds of contractions, shared
transpose DAG G⊤ L forward DAG GL compile
shared intermediates computed once
= condensed expression
subexpressions) mirrors Goldstone/MBPT diagrams
-mode
reverse
path-op
t
= relaxation: Λ / Z / ζ
optimal path DP over linear extensions; reused: grid × methods
Figure 3: The contraction-DAG leverage. The per-method term list compiles into one forward DAG that condenses shared intermediates; its transpose yields the relaxation equations for free (for EOM-CC: the ζ amplitude-response and the interstate Z-vector), and its onceoptimized path is reused across the Laplace grid and across methods. The same DAG drives the GPU code generation of Sect. 2.4. 22
Relaxation is the DAG transpose. Intuition: computing a gradient by hand means deriving a separate “relaxation” (Lagrange-multiplier) equation for every state — tedious and error-prone. This is exactly the situation automatic differentiation solves: if the energy is a computation graph (a DAG), running that graph backwards mechanically produces all the response equations, just as backpropagation produces gradients in machine learning. This view is now standard in differentiable quantum chemistry — from automatic differentiation through Hartree–Fock 10 and the PySCFAD framework 11 to the reverse-mode coupled-cluster response of Zhang et al., who obtain the relaxed (Lagrangian) response density from the backward pass instead of solving the Z-vector equation explicitly. 12 We exploit the same principle, specialized to the non-Hermitian EOM-CC response: the relaxation equations are the transpose of the density-build DAG, generated once by the machine rather than per state by hand. Standard EOM-CC analytic gradients (Stanton and Gauss; 3,4 the spin-conserving and spin-flip implementation of Krylov and co-workers 5 ) solve a left-state Λ-like Lagrange equation (with Λ the ground-state Lagrange multiplier) plus a separate orbital Z-vector (the single linear solve that yields all orbital-response contributions), each re-derived per state. Here the amplitude response is a single self-consistent multiplier ζ obtained as the transpose of the forward density-contraction DAG — identical in form across EE/IP/EA and across truncation level — and the orbital Z-vector is the same J/K build with a different density argument. Our contribution is not a new physical quantity but a mechanization: all relaxation equations are one reverse-mode pass over the forward DAG, never re-derived per state. The reverse pass is lines 11–13 of Algorithm 1. Concretely, holding L, R frozen, the amplitudes T still respond to nuclear motion; that response solves A⊤ ζ = − ξ,
ξµ =
∂F AB , ∂tµ
(22)
where A is the EOM Jacobian of Eq. (5) and ξ is the source vector ∂F/∂t, built by an analytic term-list adjoint (the reverse-mode pass) and validated against reverse-mode au23
tomatic differentiation to machine precision. It is self-consistent because ξ is the gradient of the very functional F whose densities feed the build, unifying the ground-state Λ and the excited-state response into one multiplier over the full amplitude space—replacing the per-state Z-vector prescription. A subtlety here is load-bearing for correctness: ξ is the t-derivative of the density functional F AB of Eq. (7), not of the connected matrix element ω̃(t) = ⟨L|A(t)|R⟩. The two differ P by µ (∂Ωµ /∂t)⟨L|R τ̂µ ⟩ (with Ωµ the cluster residuals; nonzero only in the single-excitation t1 /ov channel of Eq. (7), the only disconnected channel surviving at convergence), and this is precisely the disconnected t-dependence that Eq. (7) was built to carry. Using the connected form instead leaves an O(10−3 ) error in the gradient—choosing the density-functional derivative is what closes the finite-difference gate to ∼ 10−8 . Equation (22) is not a separate derivation: it is the forward density DAG run backwards (Fig. 4). What makes this mechanical, rather than a per-case insight, is that every node of GL is a single multilinear contraction, and multilinear contraction has one adjoint rule. Each forward term y += c einsum(subs; x) has the adjoint x̄ += c einsum(subs∗ ; ȳ) obtained by swapping the input and output subscript strings and accumulating into the input adjoint — the one operational rule that turns “transpose = relaxation” into a reproducible, mechanizable pass. Applied edge by edge to GL it sends every forward operation to its adjoint (J ⊤ = J, K ⊤ = K, back-transform ↔ forward-transform, projection self-adjoint) and turns the density build into the ζ solve. The result is structurally guaranteed, not fortuitous: reverse-mode differentiation of any DAG returns the transposed Jacobian-vector product of that exact DAG, so the multiplier equation it emits is necessarily A⊤ ζ = −ξ with the same A that the forward build linearizes. The same construction runs unchanged for the interstate Z-vector of the NACME (mixed L̂A , R̂B in place of L̂k , R̂k ) and for EE/IP/EA at any truncation level, because all of them are the same forward graph with different leaf tensors — one reverse pass, never a family of hand derivations. That the adjoint rule applies cleanly term by term is a deliberate consequence of how
24
the working equations are emitted. The residuals are generated symbolically (Sect. 3 below) as flat lists of single contraction terms, each carrying its coefficient, integral tags, and amplitude tags. We disable the compressed (paired) antisymmetrizers and keep only the plain transposition P (p, q), which the translator expands as 1 − swap on the external axes. Every emitted term is then one sign-explicit einsum, so the adjoint of Eqs. (22) above follows from the subscript-swap rule alone, with no antisymmetrizer sign convention to decode. Generated in this form, the same single-term structure that builds the forward density DAG and the residuals Ω̂ transposes, term for term, into the relaxation solve; emitting transpositions rather than compressed antisymmetrizers is what keeps that transpose mechanical end to end. This term-level adjoint rule also delivers the single-kernel payoff that the title asserts. Reading the rule on a J-contraction shows that a forward edge building an intermediate from a density becomes a reversed edge building the adjoint density by the same contraction; only the density argument is swapped. Every relaxation J/K build therefore shares the exact kernel signature of its forward counterpart, so one ordered AO-LT J/K kernel (Sect. 2.4) serves the forward densities and the relaxation alike — the kernel-reuse claim is a consequence of the adjoint rule, established here by construction (its measured GPU throughput is deferred to Sect. 4.4). Figure 5 is Algorithm 1 drawn as the full EOM-CCSDTQ computation graph: the forward pass (solid) is lines 5–10, the reverse-mode transpose (dashed) is lines 11–13, and both feed the single AO-LT J/K kernel that returns the gradient or NACME.
2.4
Non-symmetric J/K for transition densities
The LT-AO kernel of Sect. 2.2 is written for two distinct density arguments because the EOM transition densities demand it. We close the theory by making that requirement explicit and reducing the two-particle contraction to the one-particle J/K builds the kernel actually performs. 25
forward
J-build
t
γ
⟨f, ·⟩ + 12 ⟨g, ·⟩
F AB seed
transpose
ξ = ∂F/∂t
same J
γ̄ = ∂F/∂γ
F̄ = 1
A⊤ ζ = −ξ
ζ
Figure 4: The relaxation is automatic differentiation applied to a density build: reading the forward graph backwards produces the response equations mechanically, the same way backpropagation produces gradients in machine learning. Forward (violet): the build t → γ → F AB . Transpose (orange): its reverse-mode pass seeds F̄ = 1, propagates γ̄ = ∂F/∂γ, and reaches the source ξ = ∂F/∂t through the same J-build kernel (teal, both directions) with the adjoint density as argument; the amplitude response then solves A⊤ ζ = −ξ. Identical kernel signatures forward and backward are why one GPU kernel (Sect. 2.4) serves both. The two-particle gradient term is 12
AB x pqrs Γpqrs (pq|rs) .
P
To carry it through the one-
particle kernel we cumulant-decompose the two-particle (transition) density into its separable (mean-field) part and the connected cumulant,
c,AB AB AB 1 AB AB ΓAB pqrs = γpq γrs − 2 γps γrq + Γpqrs , {z } |
(23)
γ∧γ
written in the spatial closed-shell convention of the Notation (the 12 exchange coefficient is the closed-shell spin-traced form of the unit-coefficient spin-orbital cumulant). The separable part collapses exactly onto the ordered one-particle derivative builds of Eq. (18), with the exchange and functional factors kept explicit,
1 2
X (γ ∧ γ)pqrs (pq|rs)x = 12 J x (A, B) − 14 K x (A, B) = 12 J x (A, B) − 12 K x (A, B) ,
(24)
pqrs
where the leading 12 is the functional prefactor of 12
P
ΓAB (pq|rs)x , the inner 21 on K x is the
exchange coefficient, and J x , K x are the bare builds of Eq. (18). The two density factors
26
h, (pq|rs), hx , S x
integrals
CCSDTQ T
EOM L, R
T1
T2
R̂k = (r1 , r2 , r3 , r4 )
ε
T3
T4
L̂k = (l1 , l2 , l3 , l4 )
γ AB , ΓAB (non-symmetric)
densities
transpose transpose / kernel
AO-LT J/K kernel P x x α wα [2J −K ]
ζ = (ζ1 , . . . , ζ4 ) A⊤ ζ = −ξ
gradient dEk /dx / NACME λAB x
gradient / NACME
Figure 5: Computation graph (DAG) for the EOM-CCSDTQ excited-state gradient and NACME. Solid edges are the forward pass; orange dashed edges are the reverse-mode transpose that produces the amplitude response ζ (Fig. 4). The CCSDTQ amplitudes T1 –T4 are converged once (Dk Tk = −Rk [T ]); for each target state the biorthonormal R̂k , L̂k follow from the native two-sided non-Hermitian Davidson (Sect. 2.1). The non-symmetric transition densities γ AB , ΓAB and their transpose ζ both feed the single AO-LT J/K kernel, which assembles the gradient or NACME; every J/K build—forward or transpose, symmetric or non-symmetric—reuses the one kernel of Sect. 2.4.
27
A, B here are spin-summed total densities; the same separable build appears as 2J x − K x on the per-spatial (half) densities D(1) , Ξ of the Laplace-grid gradient Eq. (20) (the factor 2 is the closed-shell spin sum, 2J − K = 2(J − 21 K)), and as J x (γ, γ) − 12 K x (γ, γ) on the total relaxed density γ in the master gradient Eq. (8). This single chain fixes the normalization across Eqs. (23), (24), (8), and (20). The decomposition leaves only a small active-only P x cumulant remainder 12 pqrs Γc,AB pqrs (pq|rs) , so a single J/K kernel covers the separable part of every method’s two-particle term, the active-only remainder handled separately. The two one-particle factors A, B in Eq. (24) are distinct objects—in the Laplace-grid assembly they are the level-1 forward density D(1;α1 ) (bra) and the combined forward+backward ket Ξ(L;α) of Eq. (16)—and each is itself non-symmetric because L ̸= R† . Hence A ̸= B, J x (A, B) ̸= J x (B, A), and neither the bra/ket swap nor the antisymmetric part 21 (A − A⊤ ) may be dropped: a symmetric-only kernel retains 12 (A + A⊤ ) and silently discards the rest, corrupting the gradient. The symmetric ground-state J/K is the special case A = B. At the MP2 level (L = 1) the two factors coincide and Eq. (24) reduces to the symmetric J x (D, D) form of Ref. 20; only at L ≥ 2 (MP3/CCSD and above), and for every EOM transition density, is the genuinely ordered J x (D(1) , Ξ(L) ) build required—which is why a tuned, ordered kernel is the object of the GPU work below.
3
Computational Details
DAG-driven, spill-bounded GPU chunking. A naive GPU port of the non-symmetric J/K spills registers and shared memory and falls far below the roofline (the peak throughput achievable for a kernel of its arithmetic intensity)—or exceeds the 8 GB budget—on commodity hardware. The contraction DAG resolves this: its graph structure is used to split each kernel into chunks sized to the device register and shared-memory limits, so the build stays compute/bandwidth-bound rather than spill-bound. The same DAG that emits the working equations and their transpose (Sect. 2.3) thus also emits the contraction schedule
28
and its chunking—one abstraction, both payoffs. Concretely, the HGP-OS ERI and derivative recurrence is itself a DAG (a verticalrecurrence tower feeding horizontal raises), and two codegen regimes are keyed to its size. For low angular momentum (s/p/d) the DAG is unrolled —fully inlined—which is fast but explodes the live register set: a (dd|dd) quartet is 7623 nodes / 1296 outputs and spills ∼ 96 KB per thread, while an (f f f f ) quartet (∼ 78,000 nodes, 10.4 MB of source) times out ptxas (> 1500 s) and is effectively uncompilable. For f shells and above the kernel is instead looped : the ladder structure of the DAG is encoded as integer index tables walked by runtime loops, so source size is O(tables) rather than O(nodes) and the result is machine-ε equal to the unrolled path through (f f f f ). Where the unrolled register pressure still dominates, the output cone is partitioned. Output-tiling cuts the DAG into sub-DAGs capped at ∼ 2048 outputs ((dd|dd) spill 96 → 34 KB, ∼ 1.55×); frontier (antichain) chunking bounds the per-chunk live frontier but the (dd|dd) peak simultaneous-live antichain (∼ 2172) is irreducible, so the per-chunk frontier is genuinely not the peak antichain and output-tiling wins. The heaviest cc-pVTZ-class buckets exceed the 512 KB per-thread local limit; their working set is placed in an explicit HBM scratch slab indexed by a fixed thread-pool slot, with no
syncthreads and no grid synchronization, so the threads stay independent and
the kernel scales linearly (constant ∼ 220 MB scratch and ∼ 50 µs/quartet as nquartet runs 512 → 4096). NVIDIA Nsight Compute showed these f and gradient global-scratch kernels to be occupancy-starved (2–4% of SM issue), not compute-bound; a shared ∼ 3 GB scratch pool with a small block size and grid-stride slot reuse then gave 2.2–4.5×. The limiter reported by Nsight, not ptxas+wall-clock, drove the partitioning choice throughout. Each thread evaluates one contracted ERI quartet (µν|κρ) through the HGP-OS recurrence and scatters it into the Fock matrices it contributes to; in the J/K build of Eq. (17) this per-quartet scatter is
Jµν += wµνκρ (µν|κρ) Dκρ ,
Kµκ += wµνκρ (µν|κρ) Dνρ , 29
(25)
where wµνκρ is the canonical-fold multiplicity of the quartet and the accumulation is by atomicAdd into the global Fock matrices, so threads need neither cooperate nor synchronize. The thread-to-work map underlying these regimes (Fig. 6) distinguishes this kernel from prior GPU integral engines: where those map a cooperating thread block to one quartet, the present kernel maps one thread to one quartet over a grid-stride pool, externalizing each thread’s working set to an HBM global-scratch slab so the threads never cooperate or synchronize, with high-L relief from the two DAG-size-keyed codegen regimes above (unrolled with output-tiling for s/p/d, looped integer-index-table walk with external scratch for f and up) rather than from per-output kernel fission (the execution-model contrast with per-quartet GPU integral kernels is drawn in full below). The absence of
syncthreads in the looped
path is what makes the threads fully independent and the kernel scale linearly, and it is the same graph abstraction that, transposed, emits the response skeletons of Sect. 2.3.
J and K as the shared computational primitive. The Coulomb (J) and exchange (K) builds and their nuclear gradients are the universal primitive of this backend: the same kernels that assemble the Hartree–Fock reference (the SCF) also assemble the skeleton/derivative builds J x (A, B), K x (A, B) of Eq. (18) that the response gradients contract. Characterizing the DAG-generated J/K kernel on a production RHF workload therefore directly characterizes the same kernel object used for the EOM-CC response skeletons; this is what licenses the RHF and Cholesky timing and robustness data below to stand in for the response build. In particular the production launcher exposes build K screened (the symmetric L = R case, the U = V special case of Fig. 7) and the non-symmetric transitiondensity path compute densities cd nonsym (with the k G3combine nonsym ( blk/ blk f) and k G2combine nonsym combine kernels) for the transition-density (L ̸= R) exchange, P K[µν] = Qm Ũt [µ, Qm] Ṽt [ν, Qm] with distinct left/right factors U ̸= V (the symmetric build is the U = V special case)—the GPU realization of the ordered J x (A, B) ̸= J x (B, A) build derived in Sect. 2.4, here measured. We state the scope explicitly: the GPU data in this
30
(pp|pp) on the GPU: one quartet → one independent thread We map one quartet to one thread—no cooperating block, no in-block reduction. ket pairs kl
bra pairs ij
one thread
Thread T evaluates the whole (pp|pp) DAG (bottom-up): VRR-A
(e0|00)m
T VRR-C (e0|p0)m skipped (canonical fold)
HRR-AB
raise b
HRR-CD
raise d
atomicAdd → J, K
Block = a bag of independent threads. No __syncthreads, no in-block reduction; each thread scatters to J, K by atomicAdd ⇒ linear scaling—in contrast to a cooperating block that sweeps the ket list and finishes with an in-block reduction.
Figure 6: The HGP-OS recurrence kernel: one CUDA thread evaluates one quartet’s full DAG (VRR-A → VRR-C → HRR-AB → HRR-CD) over a grid-stride pool with no barriers, so threads never synchronize and the kernel scales linearly. Low angular momentum (s/p/d) is inlined into registers with output-tiling; high angular momentum (f and up) is a tablewalked loop whose per-thread working set lives in a dedicated HBM global-scratch slot. Each thread evaluates the per-quartet scatter of Eq. (25). Only the realized (bra × pivot) buckets are launched. This contrasts with the cooperating-block, ket-sweeping map of the TeraChem engine 21 (Martı́nez group), on which the AO-direct J/K formulation here builds. Measured on production RHF/CD-RHF workloads; the identical non-symmetric transition-density path (compute densities cd nonsym with the k G3combine nonsym/k G2combine nonsym combine kernels) serves the EOM-CC response (Sect. 4.4).
31
section characterize the shared J/K kernel object as exercised on production RHF and CDRHF workloads at the ∼ 730-AO scale—representative of the chromophore systems for which this backend was engineered—whereas the EOM-CC response builds invoke the identical kernel entry points (compute densities cd nonsym and its k G3combine nonsym/k G2combine nonsym combine kernels), so the kernel under test is the same object in both settings. As a kernellevel in-method characterization—distinct from a full excited-state run—we wall-time the non-symmetric transition-density gradient assembly on the same RTX 4060, driven by a surrogate (model) transition density rather than by amplitudes from a converged chromophorescale EOM solve; this is reported in Sect. 4.4 and times the assembly kernels themselves, not an end-to-end EOM-CCSD gradient (that full gradient is delivered at chromophore scale in Sect. 4.3). It anchors the characterization here at the level of the kernels actually invoked. A second design choice underlies these builds: rather than form the four-index ERI tensor per quartet, the backend Cholesky-decomposes the two-electron operator into a compact factor B and reads J and K off B as cuBLAS contraction DAGs (Fig. 7)—the N 4 tensor is never materialized. The factorization writes the ERI as
(µν|κρ) ≈
X
BQ,µν BQ,κρ ,
(26)
Q
and the Coulomb build then reads off B as two contractions over the auxiliary index Q,
Jµν =
X
BQ,µν γQ ,
γQ =
X κρ
Q
BQ,κρ Dκρ ,
h while the symmetric exchange build factors through the half-transform BQ,µi =
(27)
P
κ Cκi BQ,µκ
onto the occupied MO coefficients C (an nao ×nocc block), Kµν =
X
h h BQ,µi BQ,νi .
(28)
Q,i
For the non-symmetric (transition/interstate, L ̸= R) density the left and right half-transforms 32
differ, B U = U ·B and B V = V ·B, giving Kµν =
X
U V BQ,µm BQ,νm ,
(29)
Q,m
of which the symmetric build of Eq. (28) is the special case U = V . This compact-factorB machinery adapts the GPU J/K-over-pseudodensity construction of Hohenstein et al., 20 part of the broader low-rank two-electron-factorization program of Hohenstein, Parrish, and Martı́nez: 27 reading J and K off a low-rank factor of the two-electron operator rather than from the explicit N 4 tensor is the methodological link between that Stanford-lineage machinery and the present excited-state-response application, in which the same factored contraction is driven by left/right transition densities U ̸= V . What the present work contributes on top of that inheritance is the ordered, non-symmetric build itself: the same cuBLAS GEMM/GEMV schedule serves the symmetric SCF J/K (U = V ), the non-symmetric transition-density build J x (A, B) ̸= J x (B, A) that the non-Hermitian response densities force, and (below) the gradient—all within 8 GB.
Hardware, software, and the FP32-mixed boundary. All GPU measurements were performed on a single consumer NVIDIA GeForce RTX 4060 (Ada, sm 89, 8 GB), whose FP64 throughput is ≈ 1/64 of its FP32 rate; the J/K kernels are generated as CUDA from the contraction DAG and evaluate the AO Gaussian integrals with established GPU integral methods. 20,21 Kernels are built with -arch=sm 89 -O3 -lineinfo -use fast math; register usage and occupancy were read from ptxas -v and NVIDIA Nsight Compute. The reference (CPU) verification stack uses a spin-orbital implementation and a C backend driven from Python. Precision is mixed along a deliberate boundary. The integral-derivative and exchange builds run in FP32, while the Cholesky/CD-vector build, the two-center gradient term, the gradient densities, and the two-center metric stay FP64. The Coulomb (J) build is genuinely FP64: it is bandwidth-bound (the roofline of Sect. 4.4), so FP32 would buy nothing. “FP32-mixed” thus denotes K and the three-center gradient in FP32 with 33
One K aux-block: the N 4 quartet grid collapses into the factor B No per-quartet exchange kernel—the Cholesky factor B replaces the N 4 ERIs; per iteration K is a streamed cuBLAS GEMM over B, one aux-block Q at a time. B [ ncd ×nsig ]
(µν|λσ)
CD
host loop over Q
✗ never formed
Non-symmetric K (U ̸= V , transition density): K = B⊤ V BU with BU = UB, BV = V B. The symmetric build is the case U = V . aux-block bQ
K
One aux-block Q, four GPU steps: 1. unpack Bd [bQ , nao , nao ] cuBLAS GEMM
2. half-transform Bh = C Bd custom kernel
3. permute [Q, µ, i] → [µ, Q, i]
4. accumulate K += B⊤ h Bh (β)
Single stream, FP32 CUDA cores (no tensor cores). The cross-block sum is the cuBLAS βaccumulate, not a reduction kernel.
Figure 7: Coulomb (J) and exchange (K) energy build as cuBLAS contraction-DAGs over the compact Cholesky factor B, with the N 4 ERI tensor never formed. J is two bandwidthbound FP64 GEMVs over B; K streams the auxiliary index in ∼ 256 MB blocks as an FP32 GEMM cascade (unpack → half-transform → permute → K += B⊤h Bh ). The transitiondensity (EOM, L ̸= R) variant uses two distinct half-transforms BU = U B, BV = V B ⇒ K = B⊤V BU , of which the symmetric ground-state build is the special case U = V . The panel evaluates Eqs. (26)–(29) (J, Eq. (27); symmetric K, Eq. (28); non-symmetric K, Eq. (29)). Markers: • symmetric / ground-state density, ▲ non-symmetric / transition density; orange denotes FP64, blue FP32. Measured on production RHF/CD-RHF workloads; the identical non-symmetric transition-density path (compute densities cd nonsym with the k G3combine nonsym/k G2combine nonsym combine kernels) serves the EOM-CC response (Sect. 4.4).
34
J, the two-center gradient, and all densities in FP64. This density-precision statement holds for the symmetric ground-state path; the non-symmetric packed-response build is the one exception—there the transition-density V-transform GEMMs (the two exchange Ksandwiches and the B ·W , B⊤ ·bgK contractions) run FP32, and only the metric core (the pivots Bpiv , the pseudo-inverse J + , the half-transform W , the metric-folded Φ/Ψ, and the two-center assembly) stays FP64. A further enabling rule, given the 1/64 FP64 penalty, evaluates the Boys special function in FP32 inside an otherwise-FP64 kernel. Because the finite-difference gates of Table 1 exercise the dense FP64 reference contraction rather than this FP32-mixed GPU path, the precision boundary is validated separately: the FP32-mixed kernels reproduce their own FP64 build to machine-ε and agree with our direct four-center Cartesian ERIs to ≤ 10−7 (the same bit-exact gate detailed in the screening-robustness paragraph below), so the boundary is measured, not assumed.
FP32 is the precision floor. FP32 is also the floor : the TF32 tensor-core path sits below it and we do not use it for the amplitudes the response is built on. Running the groundstate CCSD particle-particle ladder—the O(o2 v 4 ) GEMMs whose converged amplitudes seed the EOM-CC amplitude response and the interstate Z-vector—in TF32 rather than FP32 perturbs the correlation energy of Mg-porphine (def2-SVP) by 1.5 × 10−3 Ha, at the edge of chemical accuracy, and the doubles by 7 × 10−4 in relative Frobenius norm (Section S8). The relative perturbation is ≈ 4.3 × 10−4 and, across the two systems we tested, nearly size-independent, so the absolute error tracks |Ecorr | and reaches chemical accuracy only at chromophore scale—a small-molecule proxy would hide it. The energy is a necessary-butnot-sufficient witness of amplitude quality: it is blind to the virtual-tail amplitudes TF32 truncates first and that gradients and small-gap NACMEs weight most heavily, where its data-dependent rounding could inject geometry-dependent jitter into the potential surface (argued, not scanned; Section S8). We therefore keep these amplitudes in FP32 and use no tensor cores. The temptation is real and portable: cuBLAS does not auto-dispatch the FP32
35
tensor path here, but an explicit CUTLASS tensor-op ladder GEMM runs ∼ 1.3–1.7× faster on this 4060 (∼ 10% end-to-end) and several-fold faster on datacenter parts (A100/H100), where TF32 is the first optimization a reader would reach for—which is precisely why we document that the banked seed amplitudes are where it must not be applied.
Bounded-footprint memory discipline. Device memory is claimed and released in small increments, so resident usage tracks a bounded sawtooth even for the largest systems (∼ 730 AO on the 8 GB card). This is the composition of aux-block chunking that bounds the transient dense slab to ∼ 256 MB, a growable bump-pool reused (not churned) across buckets, phase-boundary releases (the J/K buffers are freed after the SCF and before the gradient; the four-center-ERI build pool is freed once the Cholesky factor exists), JIT cubin release after the build, and never-dense invariants—only the buckets actually launched are materialized, the packed/screened layout is O(nsig naux ) rather than O(n2ao naux ), and K unpacks one aux-block at a time so the dense [naux , nao , nao ] tensor is never resident. For the 730-AO system the peak vs. resident footprint is ∼ 4.95 vs. ∼ 4.83 GB. This ceiling is for the symmetric ground-state gradient; the non-symmetric response gradient carries a higher constant—the packed transition-density relaxed density G3p [nsig , ncd ] plus a per-block scratch footprint ∼ 2× the symmetric path (its two K-sandwiches and a Φ-sandwich each need their own buffers, halving the aux-block). Even with this higher constant the responsegradient assembly is feasible for all eleven benchmark systems on the 8 GB card, because the device footprint tracks the Cholesky/integral working set—the maximum angular momentum Lmax , the Cholesky rank ncd , and the significant-pair count nsig —rather than the raw AO count. The largest system by AO count, chlorophyll a/pcseg-0 (730 AO but s, p-only, Lmax = 1, and CD-sparse), therefore carries a smaller device peak than the f-shell-dense decamethylferrocene/def2-SVP (486 AO, Lmax = 3, nsig = 163,954), the tightest case. The two largest by AO count (chlorophyll a and dna strand 4mer, 752 AO at δ = 10−3 ) fit through the device-resident G3p -eviction path that host-snapshots the FP64 B and its FP32
36
mirror to free the device for the per-pivot kernels. This assembly-feasibility question—the footprint of the gradient build—is distinct from the ground-state CCSD-solve t2 residency ceiling (Sect. 4.4; SI Section S8), which only 6 of the 11 systems clear. Both are fixed memory limits, not algorithmic ones.
Gradient execution on the same spine. The nuclear gradient runs as a two-stage DAG on the same kernel object (Fig. 8). A forward (γ → c → h → g) and an energy-weighted (Wq → U → V ) relaxed-density branch assemble the relaxed densities; being topologically independent in the DAG, the two branches run on concurrent CUDA streams. In the production FP32-mixed path the V-transform GEMMs of the transition-density build—the two exchange K-sandwiches and the B·W and B⊤·bgK contractions—run in FP32, while the metric core (the Cholesky pivots Bpiv , the pseudo-inverse J + , the half-transform W , the metricfolded Φ/Ψ, and the two-center assembly) stays FP64. Their output feeds derivative-ERI kernels generated by the L18 single-index lift (four VRR channels S0 , SA , SB , SC , with the fourth centre closed by D = −(A + B + C)); the three-center derivative kernel is FP32 while the two-center derivative and the two-center metric stay FP64. The production launcher evaluates these through a warp-per-quartet amortized 3c/2c kernel—one warp per quartet, the VRR body run lane-parallel—under
launch bounds (256,3) (3 blocks/SM, ∼ 85 reg-
isters/thread, ∼ 50% occupancy), which fills the grid that the earlier per-pivot launch left underfilled. The N = 3 bound is a sweep optimum (the three-center kernel runs 9.05 → 7.75 s, −14%; N ≥ 6 hits a register-spill cliff), giving 6.4× on the three-center build (31.2 → 4.85 s) and 3.8× on the two-center build at the f -shell scale (ferrocene/def2-SVP), bit-identical in the forces to the per-pivot path. The per-pivot ∼ 6-stream launcher is retained, now as the FP64 parity oracle: FP64 gradients run per-pivot, FP32 production runs amortized. The two-particle part of this build is the relaxed contraction X dE ∂(µν|κρ) = Γµνκρ , dX µνκρ ∂X
37
(30)
with Γ the relaxed two-particle density, assembled through the ordered J x /K x skeleton builds of Eq. (18). The centre derivatives ∂(µν|κρ)/∂X are supplied by the L18 single-index lift, which raises one bra/ket index at a time to give the three independent centre channels, the fourth closed by translational invariance, ∂ ∂ ∂ ∂ =− + + . ∂D ∂A ∂B ∂C
(31)
The gradient thus reuses the contraction-graph spine end to end: the same J x /K x builds of Eq. (18), fed the relaxed densities, with the Γ-density GEMM DAG supplying their arguments. Geometry-frozen caches for per-step reuse. The pivoted-metric pseudo-inverse J + = (B⊤piv Bpiv )+ is a function of the Cholesky pivots alone—independent of the transition densities— so it is held device-resident and reused, by a single device-to-device copy, across the many transition densities of one nuclear geometry. This skips a cuSOLVER eigensolve that is the single dominant per-call GPU cost (≈ 7.7 s/call on Mg-porphine). The cache is shared between the symmetric and non-symmetric paths (the same pivots give a bit-identical J + ) and is invalidated whenever the Cholesky factor is rebuilt—once per geometry—so it is geometryoptimization- and dynamics-safe; cache-on reproduces cache-off bit-identically, being a copy of the same buffer. A second, host-side cache reuses the geometry-frozen Schwarz/pivot screening, a marginal add-on by comparison. The combined repeated-call speedup is 1.716× (Mg-porphine) and 1.799× (dna nucleoside), almost entirely from the metric cache. This is the regime an excited-state molecular-dynamics trajectory or a multi-state response runs in—many transition densities at one fixed geometry against the same factor B—and is the cache ladder of Fig. 14(b).
Device-resident CCSD amplitudes on the same card. The response above consumes converged ground-state CC amplitudes, whose own O(o2 v 4 ) solve is the heaviest single object
38
One shell-quartet’s force: a 32-lane warp computes the L18 lift Amortized (s/p/d) path: one warp = 32 lanes per shell-quartet, cooperating with __syncwarp (never __syncthreads); lane 0 folds; atomicAdd to the forces. 1 warp = 32 lanes 0
lanes stride the VRR nodes; __syncwarp between levels (block=256, 8 warps).
forces grad[ ]
VRR: 4 derivative channels S0 , SA , SB , SC (centre 2α weights)
lane 0: 6 HRR folds raise/lower A, B, C
D = −(A + B + C) translational invariance closes the 4th centre
A
x
y
z
B
x
y
z
C
x
y
z
D
x
y
z
atomicAdd
How the warp-tasks fill the GPU: Γ-density assembly 2 independent FP64 cuBLAS-GEMM branches → 2 streams.
Per-pivot kernels 1 thread/pivot, 4–26% occupancy → roundrobin, ∼6 streams.
f -shell (Lmax ≥ 3) amortized path won’t compile → 1 thread/pivot.
Figure 8: Gradient execution. The relaxed-density assembly is a forward (γ → c → h → g) plus energy-weighted (Wq → U → V ) DAG of cuBLAS GEMMs whose two topologically independent branches run on concurrent streams; in the production FP32-mixed path the transition-density V-transform GEMMs run FP32 while the metric core (pivots, pseudoinverse, half-transform, and two-center assembly) stays FP64. Their output feeds derivativeERI kernels built by the L18 single-index lift (channels S0 , SA , SB , SC ; the fourth centre closed by D = −(A + B + C), Eq. (31)), with the three-center derivative in FP32 and the two-center derivative in FP64. Production evaluates these through a warp-per-quartet amortized 3c/2c kernel (one warp per quartet, launch bounds (256,3), ∼ 50% occupancy), 6.4×/3.8× over the earlier grid-underfilled per-pivot launch and bit-identical in the forces; the per-pivot ∼ 6-stream launcher is retained as the FP64 parity oracle. The build evaluates the relaxed two-particle contraction Eq. (30) and reuses the J x /K x kernel object of Eq. (18), so the symmetric and non-symmetric (J x (A, B) ̸= J x (B, A)) gradients differ only in the densities fed to the shared derivative kernel—the non-symmetric density assembly itself differs (the 21 (G3 + G⊤3 ) symmetrization for the canonical-pair fold and two metric-folded K-sandwiches). Measured on production RHF/CD-RHF workloads; the identical non-symmetric transition-density path (compute densities cd nonsym with the k G3combine nonsym/k G2combine nonsym combine kernels) serves the EOM-CC response (Sect. 4.4).
39
the pipeline touches. We run the closed-shell CD-RCCSD solve fully device-resident on the same 8 GB card under one organizing rule: never materialize a dense O(n4mo ) intermediate. Only the rank-3 CD factors B Q (∼ 1.6 GB) and the doubles t2 [o2 v 2 ] (∼ 3.7 GB, FP32) = stay resident; no dense oooo/ovov/oovv/ovvv/vvvv is ever formed. The residual tnew 2 D + X + X⊤ is assembled one occupied-index (i) block at a time—the full τ and the full Woooo are never built; each Woooo block is rebuilt on-device from the sliced t2,i and the occupied CD-factor blocks (Boo , Bov ). The O(o2 v 4 ) particle–particle ladder, the bottleneck, 3 is slice-streamed over the virtual index, holding only a bounded nblk a nv block in place of the
dense n4v . Pulay DIIS runs on the host—the device cannot simultaneously hold t2 , tnew 2 , and B Q at scale—so only the extrapolated amplitudes re-upload. The on-demand SCF CD engine is destroyed, freeing its ∼ 7 GB working set, before the t2 upload. On Mg-porphine (def2-SVP, 439 AO, 86 occ / 353 vir) this converges and banks the amplitudes at a 7.28 GB device peak, reproducing the SI master-table SCF/CD row bit-for-bit. The ground-state CCSD-solve wall-time campaign across the device-feasible systems—the iterative T1 /T2 cost that the Fig. 14 gradient timings exclude—is tabulated in SI Section S8, including its fp32convergence and t2 -residency limits. Fused particle–particle ladder and the device arena. A single relayout operation— the [a, e, b, f ] → [a, b, e, f ] reorder that feeds the τ -contraction—accounts for 53.6% of one ladder residual’s wall time; Nsight Compute finds its permute kernel bandwidth-bound at 201 GB/s (∼ 74% of the ∼ 272 GB/s DRAM roof) with an uncoalesced store. We delete it. A fused kernel (nla dt ladder block) folds the relayout into strided-batched cuBLAS GEMM addressing: the first GEMM, run per-a and batched over b, writes each [e, f ] block straight into its [a, b, e, f ] slot, and the second GEMM contracts τ with the transpose absorbed into the GEMM’s OP T (no separate permute kernel), emitting the [i, j, a, b] block directly (Algorithm 2). No intermediate permute kernel runs. The fused path is 1.57× faster per residual and bit-exact to the two-step path (FP64 ≤ 7 × 10−17 , FP32 ≤ 3 × 10−8 ). Every
40
transient—the G-block, the per-i slices, the einsum temporaries—is carved from one preallocated device arena by a first-fit free-list with adjacent-block coalescing; an allocation is a base-plus-offset hand-out from the slab (no driver call), not a cudaMalloc/cudaFree, so the solve neither stalls on the driver allocator nor fragments. The arena cap is sized from measured free VRAM behind a fixed margin. Consistent with the precision floor above, these amplitude GEMMs stay in FP32 and off the tensor cores. The contraction this ladder P 2 4 performs, L[ij, ab] = ef (ae|bf ) X[ij, ef ], is the dominant O(o v ) primitive not only of the ground-state residual (here X = τ2 ) but of the excited-state σ-build (with X = r2 , below): the ground-state amplitude solve is the first of two consumers of one CD-streamed, no-dense-N 4 ladder. P Algorithm 2 Fused, slice-streamed particle–particle ladder L[ij, ab] = ef (ae|bf ) τ2 [ij, ef ] on one 8 GB GPU. The [a, e, b, f ] → [a, b, e, f ] relayout is folded into strided-batched GEMM output addressing, so no permute kernel and no dense n4v tensor are ever formed. All buffers are arena allocations from the shared device pool. Require: device-resident CD factors B[Q, a, e]; τ2 [(ij), (ef )]; device arena 1: for each virtual block [a0 , a0 +na ) do ▷ slice-stream over the virtual index 2: G ← arena.alloc(na nv · nv nv ) ▷ [a, b, e, f ]; base+offset, not cudaMalloc 3: for a in the blockP do ▷ GEMM1: strided-batched over b 4: G[a, b, e, f ] ← Q B[Q, b, f ] B[Q, a, e] ▷ batched output stride writes [e, f ] into the [a, b, ·, ·] slot 5: L[(ij), (ab)] += τ2 [(ij), (ef )] G[(ab), (ef )]⊤ ▷ GEMM2: no separate permute kernel (transpose = GEMM OP T), emits [i, j, a, b] 6: arena.free(G) ▷ free-list coalesces with adjacent blocks Kernels: FP32, -arch=sm 89 (Ada, FP64 ≈ 1/64 FP32), no tensor cores
The EOM response on the same no-dense-N 4 spine. The excited-state operators inP herit the same discipline. The EOM-EE-CCSD σ = A r build routes its ef (ae|bf ) r2 [ij, ef ] ladder—over the trial-vector doubles r2 , not the ground-state τ2 —through the device CD engine (the batched non-symmetric compute densities cd nonsym path), column-chunked over the o2 trial columns at peak device memory O(chunk·n2ao ) with no dense vvvv; the ovvvintegral (singles–doubles coupling) terms stay on the host B Q stream at O(nB n2v ), no dense P ovvv. This is the same dominant O(o2 v 4 ) ladder primitive L[ij, ab] = ef (ae|bf ) X[ij, ef ] 41
as the ground-state residual of Algorithm 2, under the same no-dense-N 4 CD discipline— the same primitive, not the same kernel call: the ground-state solve drives it through
nla dt ladder block sliced over the virtual index, the EOM σ-build through compute densities cd nonsy chunked over the o2 trial columns. The two differ only in (i) the contracted object—the ground state uses X = τ2 = t2 + t1 t1 , whereas EOM uses the bare trial doubles X = r2 (since R̂ is linear, σ is linear in r2 and carries no t1 t1 term), the lighter T1 -dressed ovvv pieces staying separate in both builds—and (ii) the chunking route above. One fused CD-streamed ladder thus has two consumers: the Newton/DIIS ground-state amplitude solve and the outer non-Hermitian eigensolve. Each σ = H̄ r the EOM eigensolver evaluates is exactly this ladder-dominated H̄ r build, and the eigensolver is the outer non-Hermitian iteration layer over that inner work; this cohesion is architectural—the EOM device paths remain validated only at small scale (≤ ethene/6-31G, below), in contrast to the chromophore-scale ground-state solve. The interstate ζ/Z-vector response and the orbital (CPHF) response run AO-direct through the same engine (QouterAdjointDevice, DeviceCPHFSolver), eliminating the last dense host quartic eriao [n4ao ]; a guard aborts on any rank-4 AO allocation along the P path. The two-electron force—the relaxed contraction 2 pqrs (pq|rs)x Gpq Grs = 2J x (G, G) in the bare J x of Eq. (18), with G ≡ Gao the symmetrized relaxed AO pair density and the leading factor 2 the bra–ket pair interchange (pq|rs) = (rs|pq) (not a redefinition of J x ; the exchange channel is captured by the symmetric eigendecomposition of the symmetrized P Gao below). This is the 12 pqrs Γpqrs (pq|rs)x relaxed two-particle contraction with the separable density carried as the rank-1 outer product Gpq Grs , and it matches the device force P P 2 (pr|st)x Gao = i 2J x (WL,i , WR,i ). It is built amplitude-direct from a rank-1 (WL , WR ) eigendecomposition of the symmetrized Gao rather than a dense rank-4 AO two-particle density. Symmetrizing is exact here: the derivative integral’s pair symmetry (pq|rs)x = (rs|pq)x annihilates the antisymmetric part of the pair-density outer product, so an eigendecomposition suffices even for the non-symmetric transition densities. One ordered-argument J x kernel then serves both the symmetric ground-state (A = B) and the non-symmetric
42
(J x (A, B) ̸= J x (B, A)) builds. These response paths are free of any dense O(N 4 ) allocation by construction and are validated against the determinant oracle and finite differences at small scale (≤ ethene/6-31G; gradient and NACME against finite difference at H2 O/STO3G). At chromophore scale (Mg-porphine) these response paths are executed end to end—the complete per-atom gradient and Q–B NACME reported in Sect. 4.3 (Table 4)—but they are not validated there by an end-to-end finite-difference check, which is infeasible at that size. Throughout, we distinguish demonstrated/executed at scale from validated : end-to-end finite-difference and independent cross-code validation are established at small scale on the identical kernels (Sect. 4.3), and the chromophore-scale result is anchored only piece-wise (kernel identity, frozen-core finite differences, and machine-zero translational invariance).
Execution model: lineage and independent extensions. Figures 6, 7, and 8 together make the execution model of this backend explicit and place it within the lineage it inherits. As noted above, the AO-native, GPU-accelerated J/K-over-compact-factor machinery here inherits the J/K-pseudodensity construction of Hohenstein et al. 20 and the low-rank twoelectron-factorization program of Hohenstein, Parrish, and Martı́nez, 27 together with the GPU Gaussian-integral engineering of the TeraChem line; 21 the present application puts that inherited machinery to a response-/propagator-ready non-symmetric J/K and gradient on commodity hardware. What the present work builds independently on that foundation is fourfold. The HGPOS (Head–Gordon–Pople Obara–Saika) ERI and derivative recurrence engine is implemented independently as a Python code generator that emits CUDA kernel source, JIT-compiled at runtime through NVRTC behind a C++/CUDA host launcher, and is distinct from the McMurchie–Davidson recurrence of the TeraChem line. The ordered, non-symmetric transition-density build J x (A, B) ̸= J x (B, A)—forced by the non-Hermitian EOM/response densities and absent from symmetric-SCF J/K engines, which never need it—is realized here
as a single GPU kernel object shared with the symmetric and gradient builds (compute densities cd nonsy
43
Fig. 7); non-symmetric Cholesky-decomposed EOM-CC gradients are themselves established 28 in a conventional, MO-tensor, datacenter-memory setting; what is new here is not that gradient theory but its realization as a single AO-direct, memory-bounded GPU kernel object—the ordered build J x (A, B) ̸= J x (B, A), never forming a four-index MO tensor and shared with the symmetric SCF and gradient builds—that holds the entire non-Hermitian response within the 8 GB envelope of a consumer card, a commodity-hardware regime 28 does not target. The single contraction DAG is the natural algebraic language of the backend: one graph emits the working equations, their reverse-mode transpose, and the chunked schedule, so the symmetric SCF J/K, the non-symmetric build, and the gradient are all expressed in one abstraction (Figs. 6–8). And the whole realization is engineered for a bounded footprint on a single consumer RTX 4060 (8 GB), democratizing builds that otherwise call for datacenter accelerators. The two routes also make different, complementary design choices, each suited to its goal. The TeraChem kernels map a cooperating thread block to a quartet and close with an in-block reduction; the value kernels here map one independent thread to one quartet over a grid-stride pool and scatter straight to the global Fock matrix by atomicAdd (Fig. 6), which together with never materializing the N 4 ERI tensor—J and K are read off the compact Cholesky factor B as cuBLAS contraction-DAGs (Fig. 7)—is what keeps the build within 8 GB and linearly scaling for the commodity-hardware response workloads it targets; the TeraChem mapping is correspondingly well-matched to the high-throughput, datacenter setting it was designed for. These are different points in the design space, not a ranking. The thesis is that one abstraction—the contraction graph—serves all three of the symmetric SCF J/K (Fig. 7, U = V ), the non-symmetric transition-density build J x (A, B) ̸= J x (B, A) forced by the EOM left/right densities, and the relaxed gradient (Fig. 8)—all within an 8 GB budget.
44
Screening-collapse robustness (rediscovery of Almlöf 1982). The warning of Almlöf, Faegri, and Korsell concerns density-weighted integral screening in direct SCF 29 —which is exactly the O(N 2 ) Schwarz-screened AO-direct J/K build this paper’s kernel performs (the screened build of Sect. 2.2, whose asymptotic cost holds “once Schwarz screening is applied”). The hazard is therefore intrinsic to this work’s own evaluation path: density-weighted screening can drop Fock contributions in regions where the SCF then builds up charge, driving the energy variationally below true Hartree–Fock for any preset threshold. The same hazard reappears, in sharper form, in the Cholesky/RI factorization through which the production J/K primitive optionally runs—a loss of positive-semidefiniteness in the Choleskydecomposed (CD) two-electron operator. When the pivoted-Cholesky column-fill Schwarz screen screen tol is coupled to the decomposition threshold δ (the naive screen tol= δ), the B columns are under-filled, the CD-reconstructed operator loses positive-semidefiniteness, and RHF collapses below the ground state—by a basis-diffuseness-dependent factor (≈ 1.9× for a compact case up to ≈ 27× for a diffuse triple-ζ case). It hides because translational P invariance (| F |) and Tr[DS] = N are machine-ε on any density—exactly Almlöf’s point that only the energy exposes it. The fix is robust to arbitrary angular momentum L, for energies and forces alike. First, the screen is decoupled from the threshold, screen tol= 0.01 δ, the collapse-free knee that preserves PSD for the most diffuse manuscript basis while still fitting 8 GB (0.1 δ is insufficient for f -shell/diffuse bases and 10−10 OOMs the largest systems; only the count of significant pairs nsig grows, the Cholesky rank ncd is unchanged). Second, a theory-grounded per-step guard polices the reconstruction defect in the densityP P weighted metric, η(D) = µν ∆[µν] D[µν]2 with ∆ = Q2 − P L2 the per-pair Schwarzdiagonal completeness defect, bounded by η ≤ 21 δ Nelec (empirically η ≈ 471 collapsed vs. ≈ 2.5 × 10−4 correct against a budget ≈ 6.6 × 10−3 ). Crucially the build-time completeness certifier (defectmax ≤ δ) does not fire—only the density-weighted η does—and η is L-agnostic because ∆ is per-pair; production auto-tightens (start at 0.01 δ, rebuild ×0.01 tighter until η ≤ 21 δN ). Third, the kernels stay pure-Cartesian and the Cartesian→spherical transform is 45
applied at the matrix level (Dcart = T Dsph T⊤ in, Jsph = T⊤ Jcart T out), preserving the PSD metric for any L without touching the recurrence. In short: the density-weighted-screening collapse Almlöf warned of in direct SCF reappears as a loss of positive-semidefiniteness in the Cholesky-reconstructed operator, and decoupling the screen from δ while policing the P L-agnostic, geometry-respecting defect η = ∆D2 ≤ 21 δN restores robustness for energies and forces at arbitrary L on a single 8 GB GPU. The correctness gate is bit-exact agreeP ment with our own direct four-center Cartesian ERIs (≤ 10−7 ), with | F | ≤∼ 10−11 the translational-invariance check (necessary but not sufficient—it survived the collapse); absolute CD-RHF energies are not directly comparable to spherical references, a basis-cardinality (6d/10f ) effect rather than an error.
Basis sets and systems. The finite-difference verification uses the minimal STO-3G basis on the small closed- and open-shell systems of Table 1 (H2 O, C2 H4 , NH2 , BH2 ), spanning singlet, doublet, triplet, and quartet multiplicities.
Finite-difference protocol. Analytic gradients and interstate NACMEs are compared against central finite differences of the EOM-CCSD excitation energies, using the true biorthonormal left eigenvector Lk (not the Lk = Rk shortcut). Each Cartesian component is displaced symmetrically, and the analytic–numerical difference is gated at < 10−5 hartree/bohr.
Validation stack. Every response-density block is checked against a determinant (Fockspace) oracle and a spin-orbital reference, the ζ source vector against reverse-mode automatic differentiation, and the assembled gradient/NACME against finite differences; and against an independent code (Psi4; SI Section S5) at small and aromatic chromophore-class scale—reproducing our EOM-CCSD excitation energies to ≤ 1.5 × 10−7 Eh (to 3 × 10−8 at benzene, and to 5.6 × 10−8 code-vs-code at benzene/cc-pVDZ, 114 AO) and the excitedstate gradient to 4.6 × 10−7 Eh /a0 at H2 O and 1.4 × 10−7 at benzene—with CFOUR 28 46
the natural gold-standard extension. The oracle represents configurations as occupation bitstrings, using population-count and parity bit-operations for the Slater–Condon couplings. 30 The spin-orbital coupled-cluster amplitude residual equations through quadruples (Ω̂1 –Ω̂4 , CCSDTQ) used by the high-rank reference engine were generated symbolically with the p† q package 23,24 and consumed as differentiable einsum contractions in the single-transposition form of Sect. 2.3. The one- (Fock) and two-electron (fluctuation) parts are generated in separate passes and concatenated, which roughly halves peak memory and keeps the quadruples (Ω̂4 ) generation tractable. The harnesses are listed in Sect. 4.
4
Results and Discussion
4.1
Verification
Because the response equations are machine-generated rather than hand-derived, independent ground truth is essential; we use three. Symbolic/exact: the one-particle density γ1 and nearly all two-particle blocks reproduce a determinant (Fock-space) oracle and a spin-orbital reference bit-exactly (e.g. the chemist/physicist 2-RDM prefactor bridge to 0.0; the C backend reproduces the reference to 10−12 ). The one exception is the oovv two-particle block, where a connected-triples (T 3 ) tail leaves a ∼ 10−6 residual at the CCSD-truncation level—a truncation effect, not a coding discrepancy. Finite differences: the assembled excitationenergy gradient and interstate NACME agree with central finite differences across all four spin multiplicities (Table 1), with the true biorthonormal Lk . The ground-state CCSD gradient closes to 6.8 × 10−7 and the ζ amplitude-response solve converges by GMRES (relative tolerance 10−10 , gated at 10−7 ). These checks exercise the dense (canonical) contraction path that defines the reference densities; the convergence of the Laplace-grid realization of Sect. 2.2 with the number of quadrature points nτ is a separate axis, reported in the grid-convergence panel below.
47
Table 1: EOM-CCSD excited-state gradient and interstate NACME vs. central finite differences (STO-3G)—an implementation-correctness check (analytic assembly vs. the numerical derivative of the same code’s energy). All gated entries pass < 10−5 ; quartet NACMEs are reported best-effort (partner-availability/near-degeneracy limited). Independently reproduced via verif fegk/test multiplicities.py (all gates pass). System
multiplicity
H2 O H2 O C 2 H4 NH2 NH2 BH2 BH2
singlet triplet singlet doublet quartet doublet quartet
gradient max |A − FD|
†
1.83 × 10−6 1.44 × 10−6 8.56 × 10−7 1.97 × 10−8 5.67 × 10−7 4.08 × 10−8 2.44 × 10−8
NACME max |A − FD| 4.81 × 10−6 8.75 × 10−6 3.57 × 10−8 8.70 × 10−6 4.37 × 10−5 † 1.19 × 10−8 2.41 × 10−5 †
reported, not gated.
Comparison with exact full CI. Whereas Table 1 certifies that the analytic assembly reproduces the derivative of the energy this code computes (implementation correctness), the following test certifies that the method reproduces the exact answer (accuracy). The determinant oracle provides an exact, finite-difference-free reference for the EOMCC gradients themselves. Run at rank two it is EOM-CCSD; run at full excitation rank (K = min(Nocc , Nvirt )) it diagonalizes H̄ over the entire excitation manifold and is therefore full configuration interaction. Both gradients are assembled analytically through the same AO-LT contraction, so their difference max |gEOM-CCSD − gFCI | is the intrinsic singles-anddoubles truncation error of EOM-CCSD relative to the exact one—no finite differences enter (Table 2). For H2 the doubles space is already complete, so EOM-CCSD coincides with FCI and the analytic gradients agree identically (to < 10−12 ), a stringent exact cross-check of the gradient machinery; because both sides share the AO-LT contraction, any code-level error would surface here rather than in the truncation residuals below. For the correlated polyatomics the EOM-CCSD analytic gradients recover the exact FCI gradients to within the expected singles-and-doubles error, ∼ 10−3 Eh /a0 , with excitation energies within ∼ 10−3 Eh . The larger Table 1 system C2 H4 is beyond a dense diagonalizer—its FCI space holds 3.0×107 determinants—and remains covered by the finite-difference gates of Table 1.
48
Table 2: EOM-CCSD excited-state analytic gradients versus exact full CI (verif fegk/test fci gradient.py), with no finite differences. FockCC at rank two is EOM-CCSD and at full rank (K = min(Nocc , Nvirt )) is FCI; both gradients use the same AO-LT contraction. ∆ω is the EOM-CCSD excitation-energy error and the last column the maximum gradient deviation, both relative to exact FCI. These residuals are the intrinsic singles-and-doubles truncation error of the EOM-CCSD method, not an implementation error: because the analytic and exact gradients share the identical AO-LT contraction, any code-level error would appear in the H2 cross-check (< 10−12 ), not here. The FCI reference gradient is itself confirmed against finite differences of the FCI energy (< 10−5 ) for H2 and LiH. System H2 /6-31G LiH/STO-3G H2 O/STO-3G H2 O/STO-3G NH2 /STO-3G BH2 /STO-3G BH2 /STO-3G †
multiplicity singlet singlet singlet triplet doublet doublet quartet
FCI dim 28 495 1001 1001 2002 3432 3432
∆ω (Eh ) 0†
6.0 × 10−5 1.0 × 10−3 2.5 × 10−3 3.4 × 10−4 6.3 × 10−4 5.4 × 10−4
max |gEOM-CCSD − gFCI | < 10−12 † 1.4 × 10−4 1.8 × 10−3 1.8 × 10−2 1.9 × 10−3 2.7 × 10−3 1.5 × 10−3
EOM-CCSD ≡ FCI for two electrons (the doubles space is complete); the analytic gradients agree exactly, isolating any code-level error.
4.2
Spin-adapted open-shell EOM-CCSD: multiplicity-general gradients and NACMEs
The doublet and quartet rows of Table 1 were obtained on the spin-orbital singles-anddoubles manifold of the open-shell reference. That manifold is not closed under Ŝ 2 : 15 its excited roots carry a spin contaminant that pulls ⟨Ŝ 2 ⟩ by ∼ 3 × 10−3 from the target eigenvalue and leaves the quartet NACMEs only best-effort below the gate. Those Table 1 rows are nonetheless internally consistent—the FD gate compares the analytic assembly against the finite-difference derivative of the same fully spin-orbital state and density, so it passes regardless of spin contamination. The failure that motivates the spin-adapted path appears only when the spin-pure CSF eigenvector is paired with the singles-and-doubles-truncated density (a 1.7 × 10−3 gradient error), which the complete density below removes. We remove the contaminant at its source by solving the EOM equations in a spin-adapted, genealogically coupled configuration-state-function (CSF) basis on a semicanonical ROHF reference.
49
Each CSF is an exact eigenfunction of Ŝ 2 , so every excited root is spin-pure by construction; on the host reference path ⟨Ŝ 2 ⟩ matches the target eigenvalue to ∼ 10−15 . Spin adaptation alone does not close the gradient. The state-specific one- and two-particle densities (the unified-relaxation densities of Sect. 2.3) were originally accumulated over the singles-and-doubles excitation list, which silently discards the rank-≥ 3 determinant content that a spin-pure open-shell CSF state carries. Reinstating the complete CSF-basis state densities—the full determinant-space γ and Γ of the resolved CSF eigenvector—restores that content and closes the adiabatic-state gradient. Table 3 reports the result: spinadapted excited-state gradients and interstate NACMEs for the open-shell doublet and quartet states match central finite differences to ≤ 7.9 × 10−7 Eh /a0 (gradients) and ≤ 3.6 × 10−7 (NACMEs), passing a < 10−6 gate—including the quartet NACMEs that Table 1 could only report best-effort. The eigenvalues reproduce a full-CI-limit determinant oracle to ≤ 2.4 × 10−4 Eh . The interstate NACME numerator is built from the same complete-density construction, so the Table 3 couplings are consistent with the gradients by construction. Three reductions confirm the construction. A closed-shell singlet run, where the spinadaptation transform is the identity (U = I), is byte-identical to the closed-shell path. Driving the excitation energy to zero (ω → 0) recovers the ground-state gradient to 9.6×10−9 . Freezing the EOM eigenvectors (c, d) (holding the orbital response fixed) reproduces the Hellmann–Feynman identity to 2.41 × 10−9 , confirming that no eigenvector-response term is needed once the densities are complete. Two caveats bound the claim. First, this open-shell validation is at host gate scale (nso ≤ 18; STO-3G/6-31G)—a correctness check against finite-difference and FCI-limit oracles, not a chromophore-scale run. Second, spin purity is exact only by construction: on the device path ⟨Ŝ 2 ⟩ sits at the FP32 matvec floor (∼ 10−3 ), not the host 10−15 . The BH2 quartet gradient (7.89 × 10−7 ) sits at the finite-difference noise margin, and the degeneratemanifold gradient path (required when target states are exactly degenerate) is implemented but remains dormant and untested.
50
Table 3: Spin-adapted (CSF-basis) open-shell EOM-CCSD on a semicanonical ROHF reference. ⟨Ŝ 2 ⟩ deviation from the target eigenvalue (host reference path); excited-state gradient and interstate NACME versus central finite differences; and the excitation-energy deviation from a full-CI-limit determinant oracle. All gradients and NACMEs pass the < 10−6 Eh /a0 gate. Host gate scale (nso ≤ 18).
4.3
System
mult.
NH2 BH2 BH2
doublet doublet quartet
∆⟨Ŝ 2 ⟩
∼ 10−15 ∼ 10−15 ∼ 10−15
grad max |A − FD| 1.36 × 10−7 3.33 × 10−9 7.89 × 10−7
NACME max |A − FD| 3.6 × 10−7 1.8 × 10−8 1.2 × 10−10
|∆ω|FCI
9.7 × 10−5 1.3 × 10−4 2.4 × 10−4
Integrated device EOM-CCSD gradient and NACME on the GPU
Tables 1 and 2 validate the analytic construction on the CPU reference path. The remaining question is whether the integrated device pipeline—the GPU kernels of Sect. 3 driving the full response—reproduces that construction end to end. It does. At H2 O/STO-3G (ground-state ECCSD = −75.01228689 Eh , first excitation ω0 = 0.4572 Eh ) we ran the complete EOM-CCSD excited-state gradient and the interstate NACME numerator on the RTX 4060 through the device CD (Cholesky-decomposed) non-symmetric J/K machinery (build cd engine → device 2e force / compute densities cd nonsym → DeviceCPHFSolver) and compared against both central finite differences and the dense CPU reference (Table 6). Two scope points fix what is and is not exercised. First, at this small scale the device path uses the CD non-symmetric J/K with exact orbital-energy denominators and carries no dense N 4 AO tensor; the Laplace-grid resolution of Sect. 2.2 is not exercised here—it is validated separately as the MP2-energy denominator probe (SI Section S4) and is the chromophore-scale memory enabler. The streamed CD exchange and gradient kernels carry 4 only the Cholesky factor B Q and its pivots, with no O(NAO ) device allocation; at this small
scale, however, the host still forms a dense GAO and a dense spin-orbital response density, so the fully memory-bounded per-pair-streaming device path is a documented follow-up, not a claim of the present validation. Second, the biorthonormal EOM eigenpair (R, L) is
51
produced by the native two-sided device Davidson of Sect. 2.1, which returns both vectors from a single device-resident solve; the borrowed two-step seed (a dense Davidson right vector and a GMRES inverse-iteration left vector) is eliminated on this path and retained only as a cross-check, reproduced to ∆ω = 4.5 × 10−9 . The device-computed quantities validated here are the gradient and the NACME numerator built from this pair—the two-electron force, the generalized Fock, the ζ amplitude response (Sect. 2.3), and the CPHF orbital response. The device gradient matches central finite differences to a maximum absolute deviation of 5.05 × 10−7 Eh /a0 (1.35 × 10−6 relative) in FP64 and 1.44 × 10−6 (3.83 × 10−6 relative) in the production FP32-mixed path, both inside the < 10−5 gate; against the dense CPU reference it agrees to 4.1 × 10−14 (1.10 × 10−13 relative, i.e. FP64 round-off). The device NACME numerator ⟨LA |∂x H̄|RB ⟩ matches FD to 6.68 × 10−7 (4.02 × 10−6 relative, FP64) and 1.40 × 10−6 (8.42 × 10−6 relative, FP32), agreeing with the CPU reference to 1.9 × 10−13 (1.14 × 10−12 relative); the left/right cross pair is cleanly biorthogonal, ⟨LA |RB ⟩ = −6.2 × 10−16 . As in Sect. 2.1, the device returns only the coupling numerator; the 1/(EB − EA ) gap division and the geometric-mean symmetrization are applied afterward by the caller in an FP64 island, and same-irrep pairs (|sAB | < 10−6 ) are excluded. A control isolates the response channel: forcing the ζ amplitude-response density to zero drives the device gradient 1.55 × 10−2 Eh /a0 (4.15 × 10−2 relative) from finite differences—more than three orders of magnitude outside the gate—confirming that the ζ relaxation of Sect. 2.3 is load-bearing and that the agreement above is not accidental. This closes the validation loop at small scale: the same non-symmetric CD J/K kernel object characterized at chromophore scale in Sect. 4.4 produces, on the GPU, an EOMCCSD gradient and NACME numerator that reproduce both the finite-difference oracle and the dense reference. The integrated device EOM σ/ζ/CPHF paths are exercised end to end only at this scale (≤ ethene/6-31G); at chromophore scale the kernels and the ground-state CD-RCCSD solve are demonstrated (Sect. 3, SI Section S8), and the complete excitedstate gradient and interstate NACME are carried through to Mg-porphine below (Table 4).
52
The eigensolve itself is no longer the bottleneck at this class: the spatial-singlet production route (spin mode=’none’, Olsen preconditioning; SI Section S9) converges the six lowest singlets of benzene/STO-3G on the 4060 with the degenerate 1 E1u pair cleanly resolved (|∆ω| = 2.26 × 10−6 Eh vs. PySCF eeccsd singlet), and the converged eigenpair feeds the unchanged gradient/NACME machinery through a metric-corrected spatial→spinorbital dual embedding—an STO-3G chromophore-class demonstration. At production basis (Mg-porphine/def2-SVP; 439 AO, 86 occ / 353 vir) the full FP64 two-sided trial vector (n = 4.6 × 108 , 3.4 GiB each) does not fit alongside a biorthonormal subspace in 8 GB; a state-averaged CIS-natural-transition-orbital frozen-virtual (FNO) compression of the virtual space—built from the same Cholesky factors and the low-lying CIS amplitudes, and validated ghost-free and bit-exact against the full-space eigensolver on benzene—retains n′v = 90 of 353 virtuals and renders the FP64 eigenproblem device-resident (SI Section S10). The unchanged native two-sided Davidson, with the bright 1 Eu states selected by oscillator strength, then converges the Gouterman four-orbital manifold 31 to a drift < 5 meV: a weak Q band at 1.90 eV (f = 0.06) and a strong B/Soret band at 3.44 eV (f = 0.58), reproducing the weak-Q/strong-B intensity pattern. These bracket rather than reproduce the multireference reference: CASPT2 places Q/B at 1.78/2.65 eV 32 but is known to under estimate the porphyrin Soret, whereas EOM-CCSD in a modest basis over estimates it, so the two straddle the experimental (substituted-proxy) Soret of ∼3.0–3.2 eV from opposite sides— the large EOM-CCSD-minus-CASPT2 B gap is the expected consequence of opposite-sign method errors, not a code discrepancy. The weak Q band, where experiment, CASPT2, and EOM-CCSD all cluster near 2 eV, is the tight cross-check; only Q and B are each doubly degenerate under D4h , resolved here as clean pairs. Because this spectrum brackets rather than reproduces the multireference reference, and is obtained under an aggressive FNO truncation (Sect. 4.3), we present it as a capability demonstration—evidence that the eigensolver reaches production-basis chromophore states on the card—and not as spectroscopic accuracy evidence. No external large-scale cross-check is possible because PySCF and canonical EOM-
53
CC codes cannot run EOM-CCSD at chromophore scale on this hardware; the small-scale FD/reference/cross-code agreement (Sect. 4.3), together with the kernel-identity argument, is the validation strategy at scale.
Chromophore-scale gradient and NACME. The converged Mg-porphine eigenpair feeds the unchanged device response machinery—the mixed-frame generalized Fock, the ζ amplitude response, the CD non-symmetric two-electron force, and the closed-shell CPHF— to yield the complete per-atom EOM-CCSD excited-state gradient of the bright Q state and the Q–B interstate NACME, entirely within the 8 GB card (Table 4). Three points establish that this is a genuine end-to-end computation at chromophore scale, not a scaling extrapolation; a fourth (below) bounds the approximation it carries and is why we present it as a capability demonstration rather than a converged production result. First, the twoelectron force—the dominant cost—is the same CD non-symmetric J/K object characterized at chromophore scale in Sect. 4.4: evaluated per Cholesky-pair with a device-resident force accumulator and a J-only fast path, it completes in 17.3 min for the gradient (21,609 pairs) at a peak of 6.5 GB, directly confirming that the transition-density J/K is the critical path (density → generalized Fock 180 s → CPHF ∼62 s/coordinate → two-electron force; full perstage wall times in Table 5). Second, the frozen 29-core contribution is retained, not dropped: for the gradient it is the core-involving separable two-electron force MF(Dtot ) − MF(Dact ) (a mean-field J/K difference, 225 s) which is large—1.2 to 7.9× the active-space force, and thus essential to the total—while for the NACME the core–core block vanishes with ⟨Q|B⟩ but the core–active cross B(γ QB , Dcore ), linear in the transition density, survives and is included. Third, because no finite-difference oracle is affordable at this scale, correctness is anchored piece-wise: the integrated device gradient and NACME numerator against H2 O FD (Table 6; 9.3 × 10−12 and 4.0 × 10−6 ), and the frozen-core gradient and NACME-cross terms against their own FD ground truths (4.1 × 10−9 and 6.7 × 10−10 ); the assembled two-electron P force is then translationally invariant to machine zero (| A FA | ≤ 7 × 10−12 Eh /a0 over
54
all 37 atoms), the physical confirmation that the force is complete. The dominant uncertainty in both the chromophore-scale gradient and NACME is the FNO virtual truncation itself: a direct gradient-convergence study on affordable proxies (SI Section S10, its gradientconvergence table) finds the FNO excited-state gradient to be a convergent but non-smooth approximation—bit-exact at full retention, yet at the 90/353 ≈ 25% retention used here preserved only to ∼10−2 Eh /a0 per component, and more sensitive than the excitation energy because it weights the virtual tail the compression discards first. The NACME, built from the same compressed virtual space, inherits a truncation uncertainty of this same ∼10−2 class—not separately convergence-tested, as the study probed gradients—beneath which an additional on-device numerical floor is set by the FP32 FNO-Davidson eigenvectors, whose Q/B biorthogonality ⟨LQ |RB ⟩ = −4.6 × 10−4 (versus 10−12 at all-electron H2 O) contributes a further ∼10−3 relative. We therefore report the Mg-porphine gradient and NACME as a complete, executed, memory-bounded chromophore-scale capability demonstration carrying a stated ∼10−2 Eh /a0 FNO-truncation uncertainty—not a converged production number— and recommend a per-system retention check (recompute at a larger n′v , for both the gradient and the NACME) before any production use.
Validation strategy at chromophore scale. Two properties of this at-scale result must P be read precisely. First, the machine-zero translational invariance | A FA | ≤ 7×10−12 Eh /a0 certifies that the assembled force is complete—that every term contributing to dE/dx is summed with the correct sign and prefactor—but it does not certify that the underlying densities are correct; it is a necessary-but-not-sufficient check. Second, with no per-atom finite-difference oracle and no cross-code reference affordable at this scale, correctness rests on an explicit chain of small-scale certifications transported to Mg-porphine by kernel identity: (i) the integrated device gradient and NACME numerator reproduce H2 O finite differences at the 10−6 –10−12 level (Table 6); (ii) the frozen-core two-electron force and the NACME core–active cross are each finite-difference validated in isolation (4.1 × 10−9 , 6.7 × 10−10 );
55
(iii) the biorthogonality bound ⟨LQ |RB ⟩ = −4.6 × 10−4 caps the residual NACME preciP sion at ∼10−3 ; and (iv) the machine-zero A FA confirms force completeness. Because the chromophore-scale kernels are the same objects exercised end to end at small scale (Table 6), this chain transports the small-scale agreement to Mg-porphine. A component-wise comparison against an independent implementation (CFOUR-class) on a small molecule is the committed external cross-check and is left as future work (M4). Table 4: Complete device (RTX 4060) EOM-CCSD excited-state gradient and Q–B interstate NACME at Mg-porphine (def2-SVP, 439 AO, 29 frozen core; from the FNO-compressed EOM eigenpair, entirely within 8 GB, peak GPU ≤ 6.5 GB / host ≤ 54 GB). Bright-Q gradient deFULL = deELEC + Fact 2 + Ffc (one-electron + orbital-relaxation, active-space and frozencore two-electron force) in Eh /a0 ; NACME dQB in a−1 (Q–B gap ωB − ωQ = 1.546 eV). 0 Representative atoms shown; the full two-electron force is translationally invariant over all 37 atoms to ≤ 7 × 10−12 Eh /a0 . atom
gradient deFULL (Eh /a0 )
NACME dQB (a−1 0 )
Mg C N
(−1.381, 0.228, 0.675) (−6.351, −7.741, −1.725) (−0.233, 2.574, 0.613)
(−0.451, 0.087, 0.011) (−0.890, 2.816, −0.579) (−3.467, 4.943, 0.322)
Table 5: Per-stage wall time and peak GPU memory for the complete device EOM-CCSD excited-state gradient (bright Q) and Q–B NACME at Mg-porphine (def2-SVP, 439 AO, FNO n′v = 90) on the RTX 4060. The FNO density and the CD non-symmetric transitiondensity J/K two-electron force dominate; the two-electron force is the single largest device stage, confirming it is the critical path. Peak memory stays within the 8 GB card (host peak ≤ 54 GB); CPHF is quoted per nuclear coordinate. stage
gradient
FNO density (Γ, ξ, ζ, 2-RDM) device RHF (439 AO) generalized Fock (mixed frame) CPHF orbital response active-space 2e-force (21,609 pairs) frozen-core 2e-force
32 min
NACME
24 min 50 s 3.0 min 62 s/coord 91 s/coord 17.3 min 17.2 min 3.8 min 4.0 min
56
peak GPU 7.5 GB 2.9 GB 2.9 GB 4.1 GB 6.5 GB 4.2 GB
Table 6: Integrated device (RTX 4060) EOM-CCSD excited-state gradient and interstate NACME numerator at H2 O/STO-3G (ECCSD = −75.01228689 Eh , ω0 = 0.4572 Eh ), seeded by the native two-sided non-Hermitian Davidson and run through the device CD non-symmetric J/K path with exact denominators, versus central finite differences and the dense CPU reference. Maximum absolute deviation max |·| in Eh /a0 with the relative deviation in parentheses; the gate is < 10−5 . The ζ=0 row is a control with the amplitude-response density forced to zero. Single self-consistent run; reproduced via tests/validate eom gpu grad.py (gradient) and tests/validate eom gpu nacme.py (NACME). Device quantity
vs. FD (FP64)
vs. FD (FP32)
vs. CPU ref.
EOM-CCSD gradient
5.05 × 10−7 (1.35 × 10−6 ) 6.68 × 10−7 (4.02 × 10−6 )
1.44 × 10−6 (3.83 × 10−6 ) 1.40 × 10−6 (8.42 × 10−6 )
4.13 × 10−14 (1.10 × 10−13 ) 1.90 × 10−13 (1.14 × 10−12 )
—
—
NACME numerator gradient, ζ=0 (control)
4.4
1.55 × 10−2
Performance on the RTX 4060
We characterize the DAG-generated J/K and gradient kernels on the RTX 4060 at the scale the backend was engineered for—production Hartree–Fock and CD-RHF runs on the chromophore systems of Table 1 up to 730 Cartesian AO. Because the same kernel object serves the SCF reference and the response skeletons (Sect. 3), this characterization is of the very kernels the EOM-CC gradient invokes; the symmetric builds are the ground-state-density case and the non-symmetric build is the transition-density case (J x (A, B) ̸= J x (B, A), Sect. 2.4). The in-method response-gradient build—the wall-time of the non-symmetric transition-density gradient, its per-step (repeated-call) cost under nuclear dynamics, and its parity with the ground-state gradient—is given in Fig. 14; because the contraction cost is independent of the density values, a surrogate transition density reproduces the wall-time of a converged EOM-CCSD eigenvector. The roofline (Fig. 9) places the symmetric and non-symmetric J/K kernels and the gradient builds against the device ceilings; the per-step and per-iteration decompositions (Figs. 10–11), the scaling exponents (Fig. 12), and the speedups over single-thread PySCF (Fig. 13) characterize the shared kernel across the chromophore benchmark set. The headline 57
OSX-GPU kernel roofline RTX 4060 (AD107, sm_89) Achieved performance (GFLOP/s)
FP32 peak 15.1 TFLOP/s [FP32] K exchange (non-symmetric) 4414 GFLOP/s · 67% of bandwidth roof
104
[FP32] K exchange (symmetric) 2007 GFLOP/s · 80% of bandwidth roof
103
[FP64] J Coulomb (symmetric) 63 GFLOP/s · 94% of bandwidth roof
102 AM DR
f roo
/s GB 2 27
[FP32] CD -density (build) 212 GFLOP/s · 9% of bandwidth roof
[FP32] CD 3c-grad (sym=nonsym) 95 GFLOP/s · 5% of bandwidth roof
FP64 peak 0.236 TFLOP/s
[FP64] CD -density (build) 194 GFLOP/s · 82% of FP64 roof
[FP64] CD 2c-grad (sym=nonsym)
101 34 GFLOP/s · 14% of FP64 roof 10 1
[FP32] CD -density (response (non-symmetric)) 3679 GFLOP/s · 52% of bandwidth roof
FP32 ridge AI = 56
101
100
Arithmetic intensity (FLOP / byte)
symmetric / shared (ground-state density) non-symmetric (transition density)
102
FP32 compute / point FP64 compute / point
Figure 9: Kernel roofline on the RTX 4060 (AD107, sm 89): achieved throughput vs. arithmetic intensity for the production J/K and gradient kernels, against the FP32 (15.1 TFLOP/s) and FP64 (0.236 TFLOP/s) compute peaks and the 272 GB/s DRAM roof (FP32 ridge at AI = 56). Circles are the symmetric/shared (ground-state-density) builds; triangles are the non-symmetric transition-density builds—the GPU realization of the ordered J x (A, B) ̸= J x (B, A) build of Sect. 2.4. The Coulomb build [FP64] (63 GFLOP/s, 94% of the bandwidth roof) is bandwidth-bound, which is why it stays FP64; symmetric K [FP32] (2007 GFLOP/s, 80%) and the non-symmetric transition-density K [FP32] (4414 GFLOP/s, 67%) are likewise bandwidth-bound; the three-center gradient [FP32] (95 GFLOP/s, 5%) is latency/occupancy-bound, and the two-center gradient [FP64] (34 GFLOP/s, 14% of the FP64 roof) is the only compute-leaning kernel. The transition-density Γ-density build [FP32] (3679 GFLOP/s, 52% of the bandwidth roof at AI = 26)—the response gradient’s distinctive kernel—carries ∼ 3× the arithmetic intensity of the symmetric Γ-density build (its two K-sandwiches vs. one), lifting SM utilization from 15% to 58% and placing it well up the FP32 line; the metric pseudo-inverse is excluded as it is cached once per nuclear geometry.
58
OSX-GPU AIMD-step cost decomposition (headline chromophores) (one-time) chlorophyllCDa build · pcseg-0 (137 atoms, 730 AO)
(a)
One AIMD step (s)
40.0%
0
9.9%
SCF energy
mg porphine · def2-SVP (37 atoms, 439 AO)
(b)
50.1%
25 50 75 Share of one AIMD step (%)
Analytic gradient
113.4 s
100
47.3%
0
(c)
49.4%
25 50 75 Share of one AIMD step (%)
-density 48.6%
56.8 s
Analytic-gradient breakdown Per-SCF-iteration energy breakdown
100
(d) -density 70.8%
(e)
43.8 s
-density (70.8%) 3c-grad (20.1%)
21.6 s
2c-grad (3.0%) 1e-grad (6.2%)
(f)
-density (48.6%) 3c-grad (42.1%)
Exchange (K) 59.7%
450.4 ms
Coulomb (J) (5.6%) Exchange (K) (59.7%)
2c-grad (7.4%) 1e-grad (1.9%)
Exchange (K) 42.5%
97.6 ms
Diag (eigh) (17.8%) Rest (16.9%)
Coulomb (J) (14.0%) Exchange (K) (42.5%)
Diag (eigh) (25.6%) Rest (17.8%)
Figure 10: Cost decomposition of one Born–Oppenheimer molecular-dynamics step for the two headline chromophores, chlorophyll a (pcseg-0, 137 atoms, 730 AO) and Mg-porphine (def2-SVP, 37 atoms, 439 AO). (a,b) the one-time Cholesky build vs. the per-step analytic gradient (the gradient is 80.4% / 73.9% of the step). (c,d) the analytic gradient is dominated by the Γ-density contraction (91.8% / 80.7%), the three- and two-center gradient builds and the one-electron term making up the rest. (e,f) the per-SCF-iteration energy build is exchange-dominated (K 59.7% / 42.5%). These are the same J/K and gradient kernels the EOM-CC response invokes (Sect. 3).
59
OSX-GPU per-SCF-iteration energy-term breakdown (systems ordered by AO count) Coulomb (J)
Share of one SCF iteration (%)
(a)
100
Diag (eigh)
Rest
80 60 40 20 0
(b) Total time per SCF iteration (ms, log)
Exchange (K)
Share of one SCF iteration (%)
Total walltime per SCF iteration (ms) 102
Basis set (bar colour) def2-SVP def2-TZVP 24
23
48
62
pcseg-0
98
63
112
450 107
414
139
101 100
ne oce ) ferr236 AO (
e mer ) side ) zen ) ben52 AO nucle3o5 AO trand 726 AO 2 3 3 s a ( ( dn ( dna
e in ne obil ) phin ) oct6a AO) po4r39 AO ocya4n63 AO g 9 3 c m ( ( ( phy
r a ne ne 4me ) hyll ) oce ) oce ) ferr472 AO dmferr486 AO hloro7p30 AO strand752 AO c ( ( ( ( dna
Figure 11: Per-SCF-iteration energy-term breakdown across the eleven benchmark systems, ordered by Cartesian AO count. (a) relative shares of Coulomb (J), exchange (K), diagonalization (eigh), and the remaining work; (b) absolute wall-time per SCF iteration (log scale; bar color denotes the basis set), from 24 ms (ferrocene, 236 AO) to 450 ms (chlorophyll a, 730 AO).
60
OSX-GPU scaling on a single RTX 4060 (gradient systems, 31 (a) Per-iteration J / K and gradient scaling vs AO count
1.0
102
Cholesky sparsity vs system size
benzene
0.9
octane
ferrocene (def2-SVP)
ferrocene (def2-TZVP)
0.8 101
(b)
Gradient walltime (s)
K
(c)
N 3.1
102
3 × 102
4 × 102
6 × 102
-density grad N 2.9 total gradient N 1.4 total gradient (8 GB ceiling)
2 CD sparsity nsig/nAO
Per-iteration walltime (ms)
J
N 1.9
137 atoms)
dmferrocene
0.7 0.6 0.5
dna nucleoside
dna strand 2mer
mg porphine dna strand 4mer
0.4 101 100
open grey : 8 GB feasibility ceiling (per-pivot fallback)
3 × 102
4 × 102
Cartesian AO count N
6 × 102
0.3 0.2
other (conjugated) (saturated) 3 × 102
phycocyanobilin chlorophyll a 4 × 102 6 × 102 Cartesian AO count N
Figure 12: Scaling on a single RTX 4060 over the gradient systems (31–137 atoms). (a) periteration J and K wall-time vs. Cartesian AO count N , with fitted exponents (J ∝ N 1.9 , K ∝ N 3.1 ); (b) gradient wall-time for the eight systems using the optimized path, the device-resident FP32-mixed Γ-density build ∝ N 2.9 and the total amortized gradient ∝ N 1.4 (the warp-per-quartet 3c/2c kernels now apply at all L, including the f-shell, so the total is no longer dominated by the recurrence tower); the three largest systems at the 8 GB feasibility ceiling (per-pivot fallback) are shown separately as open markers; (c) the Cholesky sparsity nsig /n2ao vs. system size, contrasting the denser π-conjugated systems with the sparser saturated ones.
61
(a) J+K Fock build (FP32-mixed CD vs FP64 direct-4c)
(b) Analytic 2e gradient (FP32-mixed CD vs FP64 direct-4c, OMP=MKL=1)
parity (1×) ) P) P) -0) P) P) -0) P) P) -0) -0) P (SV ZV (SV eg ZV (SV eg ZV (SV eg eg Fc ne (T side (pcs ne (T orph (pcs Fc (T mFc (pcs (pcs d hl-a er ze leo er cta -p PCB C a-4m ben nuc a-2m o Mg dn dn
101
33×
33×
34× 3.0×
15×
11×
15× 1.3×
100
6.1×
43×
101
6.8×
1817×
728×
1260×
1074×
797×
644×
1296×
1182×
641×
102
550×
103
Speedup (×)
s/p/d shells (VRR-amortised) f-shell (recurrence-bound, unoptimised) parity (1×)
63×
Speedup (×)
104
P) P) P) -0) P) P) -0) P) P) -0) -0) (SV ZV (SV eg ZV (SV eg ZV (SV eg eg Fc ne (T side (pcs ne (T orph (pcs Fc (T mFc (pcs (pcs d hl-a er ze leo er cta -p PCB C a-4m ben nuc a-2m o Mg dn dn
Figure 13: Speedup of the GPU backend (OSX-GPU, CD, FP32-mixed) over single-thread PySCF (FP64 direct ERIs, OMP NUM THREADS = MKL NUM THREADS = 1), the same physical observable on both sides, across the eleven systems. This single-thread, FP64, directERI baseline confounds three effects—hardware, density fitting (Cholesky), and FP32mixed precision—so panel (a) is a reference-implementation comparison rather than a likefor-like kernel benchmark; the defensible same-machine figure is the analytic gradient of panel (b). (a) the per-SCF-iteration J+K Fock build, 63×–1817× (reference-implementation comparison); (b) the analytic two-electron gradient, 6×–43× across the eight systems that fit the 8 GB device with the optimized path (device-resident FP32-mixed Γ-density build + VRR-amortized warp-per-quartet 3c/2c kernels at all L, including the f-shell). The three feasibility-ceiling systems (open markers; two f-shell—ferrocene/def2-TZVP and decamethylferrocene/def2-SVP—plus one large pcseg-0 system, dna strand 4mer) exceed device memory with both optimizations co-resident and retain the per-pivot gradient (1.3×– 33×); this is a fixed 8 GB-memory limit, not an algorithmic one. Dashed line: parity.
62
In-method transition-density response gradient on a single RTX 4060 Per-system cost decomposition
25 20
-density (response) 3c-grad 2c-grad
Repeated-call cache ladder
(b)
off scale (device-resident path): dmferrocene/def2-SVP: 111 s
20
Per-call wall-time (s)
Response-gradient wall-time (s)
(a)
15 10 5 0
er er ne ine ne ide ne ne ilin oc2e-SVPbenz-eTZVPcleos2-SVPnd 2mseg-0 oct-aTZVPorph2-SVPanobseg-0rroc-eTZVPnd 4mseg-0 r r fe def ef2 nu def tra pc ef2 g p def ocy pc fe ef2 tra pc d dna d m hyc d as as p dn dn
no cache + metric cache + host-screen cache
1.7×
15 10
1.8×
5 0
dna nucleoside def2-SVP
mg porphine def2-SVP
Figure 14: In-method transition-density response gradient on the RTX 4060 (the production FP32-mixed CD path, the same 3c/2c kernels as the ground-state gradient). All eleven response-gradient assemblies are feasible on the 8 GB device; the two largest by AO count (chlorophyll a, 730 AO; dna strand 4mer, 752 AO) and the f-shell decamethylferrocene/def2SVP (486 AO, the tightest case, not the largest) require the device-resident eviction path. Feasibility tracks the CD footprint (Lmax , ncd , nsig ), not the AO count, and is distinct from the ground-state CCSD-solve t2 ceiling (6 of 11 systems; SI S8) that this figure excludes by construction. (a) Per-system cost: the response gradient is three-center-dominated, and the non-symmetric transition-density Γ build—the only kernel specific to the ordered J x (A, B) ̸= J x (B, A) density—is a thin slice (0.1–1.7 s). Decamethylferrocene fits only via the eviction path; its fused 3c+2c launch runs at ∼ 3.6× the per-call cost of the comparably sized ferrocene/def2-TZVP (∼ 110 s, off-scale), and its degenerate per-term split (ratio 0.757) is the price of the sub-1 GB margin, not comparable term-by-term to the standard path. (b) Repeated-call cache ladder: reusing the geometry-frozen Cholesky metric across the many transition densities of one geometry gives a 1.7–1.8× per-step speedup—the regime of an excited-state AIMD trajectory. On the standard (non-evicting) path the response gradient matches the ground-state gradient to within 1% (ratio 0.99–1.00), since one kernel serves both density symmetries and the contraction cost is density-independent. This is a cost/kernel-identity assembly timing (a surrogate transition density; resp/ground ≈ 1.00 by construction), not a chromophore-scale EOM-gradient demonstration (the full end-to-end run is M3); it complements the small-scale end-to-end device validation of Sect. 4.3.
63
same-machine figure is the analytic two-electron gradient, 6×–43× (Fig. 13(b)); the larger Fock-build ratios (63×–1817×, Fig. 13(a)) are taken against a single-thread, FP64, directERI reference and confound hardware, density fitting, and precision, so we report them as a reference-implementation comparison, not a headline. The in-method response gradient closes this loop at the kernel level (Fig. 14): the non-symmetric transition-density gradient assembly (driven by a surrogate transition density, not a chromophore-scale EOM solve, and complementing the small-scale end-to-end device validation of Sect. 4.3), wall-timed on the 4060 across all eleven feasible response-gradient assemblies, is three-center-dominated with the transition-density Γ-density on its critical path, matches the ground-state gradient to within 1% (0.99–1.00×, one kernel for both density symmetries), and inherits the 6×–43× kernel speedup of Fig. 13, with a further 1.7–1.8× from reusing the geometry-frozen Cholesky metric (and, marginally, the screening) across the per-MD-step transition densities. The Laplace-grid convergence of the AO-LT kernel—gradient/energy error versus the number of Laplace points nτ , exercising the Laplace path that the dense reference checks of Sect. 4.1 do not—is reported in SI Section S4 (Fig. S1). That study probes the shared minimax-Laplace grid through the LT-MP2 energy as a proxy on the same grid, rather than as a direct EOM-gradient-versus-nτ sweep. The DAG-derived chunking’s effect—shifting the kernel from spill-bound to compute/bandwidth-bound—is evidenced by the register-pressure and Nsight-Compute occupancy characterization above; a dedicated chunked-versus-unchunked ablation table (ptxas -v spills, occupancy, throughput) would isolate it causally and is left as a future refinement (SI Section S11).
5
Conclusions
We have realized the complete non-Hermitian EOM-CCSD excited-state gradient and interstate NACME within the 8 GB of a consumer GPU: validated end to end at small scale— including an independent cross-code (Psi4) check that reaches the excited-state gradient
64
itself and holds from H2 O to aromatic benzene—and executed at the chlorophyll-core chromophore Mg-porphine as a bounded capability demonstration. The enabling mechanism is a single contraction DAG whose reverse-mode transpose is the EOM-CC relaxation, from which the symmetric, non-symmetric, and gradient J/K builds are all emitted. By construction the DAG spine and LT-AO kernel are designed to apply across MP2– MP6, CCSD–CCSDTQ, 33 and the CASSCF/CP-MCSCF/DMRG-SCF/CASPT2 families— each filling the same S1–S5 template, differing only in which densities are built and how they are relaxed. This generality is a structural property of the shared LT-AO kernel of Sect. 2.2 and the DAG transpose of Sect. 2.3; in the present work it is the EOM-CC case that we develop and validate in full, and we do not claim the other families as tested here. In the multireference regime in particular, the cumulant collapse of Eq. (24) leaves a large active-only remainder that the single J/K kernel does not absorb, so the CASSCF/CASPT2 reach is a structural property of the template rather than a demonstrated capability. The reason this generality is structural rather than coincidental is that the construction operates on the contraction graph, which mirrors the Goldstone/MBPT diagrammatic expansion (Sect. 2.3): for a property expressible as a finite, truncated sum of diagrams and cast as a differentiable (Lagrangian) functional, the reverse-mode transpose of the diagram-contraction graph is the corresponding relaxation/Lagrange-multiplier equation, generated rather than re-derived. This is a design property of the framework, demonstrated here only for nonHermitian EOM-CC. It is not a statement about infinite resummations (e.g. RPA or ladder series) or frequency-dependent self-energies, where the response follows from differentiating the underlying self-consistency or linear-response solve rather than transposing a finite term list; extending the present mechanism to Green’s-function/self-energy formulations is therefore an open question rather than a corollary. The reach of the method is bounded by the assumptions of its two pillars. The shifted Laplace quadrature of Eq. (21) requires ω < mink Dk , restricting the excited states to the HOMO–LUMO-gap regime; states above the first denominator pole fall outside the present
65
grid. The AO-direct formulation eliminates one power of the virtual-space dimension at every truncation (the entire vvvv block at the doubles level), which—together with the spillbounded chunking—is what makes correlated excited-state gradients and NACMEs fit and run on commodity hardware. A natural extension, squarely in the same machinery, is analytic gradients/NACMEs for ab-initio cavity-QED (QED-EOM-CC) polaritonic states, where the non-symmetric transition densities and the DAG-transpose response carry over unchanged.
Data and code availability The validation stack is publicly deposited so that the paper’s accuracy claims can be independently checked against the same references: the verification suites—the determinant (Fockspace) oracle, the finite-difference gradient/NACME harness (test multiplicities.py), and the full-CI gradient harness (test fci gradient.py)—together with the independent cross-code comparison scripts (Section S5), the molecular geometries of Tables 1–2, the complete working-equation derivations, and portions of the PTX code-generator output (the emitted GPU-kernel assembly for the symmetric, non-symmetric, and gradient builds), are deposited on Zenodo (DOI assigned on acceptance) and mirrored in a public GitHub repository. These host and independent-code checks reproduce the small-scale accuracy references without the proprietary kernels, and Section 3 describes the device algorithm—the contraction schedule and its live-set chunking, the shared-memory lift, the arena allocator, and the fused particle–particle ladder—in sufficient detail to reimplement the kernels from the paper. The performance-tuned production CUDA J/K kernels—and hence the timing, roofline, and chromophore-scale execution data—are proprietary to the funder and available from the corresponding author upon reasonable request.
66
Supporting Information Available The Supporting Information is available free of charge at http://pubs.acs.org. It contains: the computational environment and hardware/software for the CPU verification stack and the GPU backend (Section S1); full Cartesian molecular geometries for all finite-difference and full-CI systems (S2); the coupled-cluster, EOM, finite-difference, and full-CI numerical protocols and tolerances (S3); the Laplace-grid (nτ ) convergence study (S4, Figure S1); expanded verification tables (S5); the determinant-oracle and p† q amplitude-generation details (S6); the complete working-equation derivations (S7); GPU-backend reproducibility, timing, the device-resident CCSD-solve campaign, and roofline data (S8); an algorithmic characterization of the non-Hermitian EOM Davidson eigensolver on synthetic small matrices (S9); the FNO virtual compression for the production-basis eigensolve (S10); the status of in-method deliverables (S11); and the data- and code-availability statement (S12).
Acknowledgement We acknowledge financial support and computational resources provided by NeuroTechNet S.A.S. The validation harnesses, cross-code comparison scripts, and molecular geometries supporting the findings of this study are publicly deposited (Zenodo, DOI on acceptance; public GitHub mirror); the production CUDA kernels are available from the corresponding author upon reasonable request. The author declares no competing financial interest.
References (1) Stanton, J. F.; Bartlett, R. J. The equation of motion coupled-cluster method. A systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited state properties. J. Chem. Phys. 1993, 98, 7029–7039.
67
(2) Shavitt, I.; Bartlett, R. J. Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory; Cambridge University Press: Cambridge, 2009. (3) Stanton, J. F. Many-body methods for excited state potential energy surfaces. I. General theory of energy gradients for the equation-of-motion coupled-cluster method. J. Chem. Phys. 1993, 99, 8840–8847. (4) Stanton, J. F.; Gauss, J. Analytic energy derivatives for the equation-of-motion coupledcluster method: Algebraic expressions, implementation, and application to the S1 state of HFCO. Theor. Chim. Acta 1995, 91, 267–289. (5) Levchenko, S. V.; Wang, T.; Krylov, A. I. Analytic gradients for the spin-conserving and spin-flipping equation-of-motion coupled-cluster models with single and double substitutions. J. Chem. Phys. 2005, 122, 224106. (6) Ichino, T.; Gauss, J.; Stanton, J. F. Quasidiabatic states described by coupled-cluster theory. J. Chem. Phys. 2009, 130, 174105. (7) Tajti, A.; Szalay, P. G. Analytic evaluation of the nonadiabatic coupling vector between excited states using equation-of-motion coupled-cluster theory. J. Chem. Phys. 2009, 131, 124104. (8) Faraji, S.; Matsika, S.; Krylov, A. I. Calculations of non-adiabatic couplings within equation-of-motion coupled-cluster framework: Theory, implementation, and validation against multi-reference methods. J. Chem. Phys. 2018, 148, 044103. (9) Kjønstad, E. F.; Koch, H. Communication: Non-adiabatic derivative coupling elements for the coupled cluster singles and doubles model. J. Chem. Phys. 2023, 158, 161106. (10) Tamayo-Mendoza, T.; Kreisbeck, C.; Lindh, R.; Aspuru-Guzik, A. Automatic differentiation in quantum chemistry with applications to fully variational Hartree–Fock. ACS Cent. Sci. 2018, 4, 559–566. 68
(11) Zhang, X.; Chan, G. K.-L. Differentiable quantum chemistry with PySCF for molecules and materials at the mean-field level and beyond. J. Chem. Phys. 2022, 157, 204801. (12) Zhang, X.; Li, C.; Ye, H.-Z.; Berkelbach, T. C.; Chan, G. K.-L. Performant automatic differentiation of local coupled cluster theories: Response properties and ab initio molecular dynamics. J. Chem. Phys. 2024, 161, 014109. (13) Hirata, S. Tensor contraction engine: Abstraction and automated parallel implementation of configuration-interaction, coupled-cluster, and many-body perturbation theories. J. Phys. Chem. A 2003, 107, 9887–9897. (14) Schweizer, S.; Doser, B.; Ochsenfeld, C. An atomic orbital-based reformulation of energy gradients in second-order Møller–Plesset perturbation theory. J. Chem. Phys. 2008, 128, 154101. (15) Pauncz, R. Spin Eigenfunctions: Construction and Use; Plenum Press: New York, 1979. (16) Almlöf, J. Elimination of energy denominators in Møller–Plesset perturbation theory by a Laplace transform approach. Chem. Phys. Lett. 1991, 181, 319–320. (17) Häser, M.; Almlöf, J. Laplace transform techniques in Møller–Plesset perturbation theory. J. Chem. Phys. 1992, 96, 489–494. (18) Takatsuka, A.; Ten-no, S.; Hackbusch, W. Minimax approximation for the decomposition of energy denominators in Laplace-transformed Møller–Plesset perturbation theories. J. Chem. Phys. 2008, 129, 044112. (19) Hackbusch, W. Computation of best L∞ exponential sums for 1/x by Remez’ algorithm. Comput. Vis. Sci. 2019, 20, 1–11. (20) Hohenstein, E. G.; Bouduban, M. E. F.; Song, C.; Luehr, N.; Ufimtsev, I. S.;
69
Martı́nez, T. J. Analytic first derivatives of complete active space self-consistent field energies on graphical processing units. J. Chem. Phys. 2015, 142, 224103. (21) Wang, Y.; Hait, D.; Johnson, K. G.; Fajen, O. J.; Zhang, J. H.; Guerrero, R. D.; Martı́nez, T. J. Extending GPU-accelerated Gaussian integrals in the TeraChem software package to f type orbitals: Implementation and applications. J. Chem. Phys. 2024, 161, 174118. (22) Aho, A. V.; Lam, M. S.; Sethi, R.; Ullman, J. D. Compilers: Principles, Techniques, and Tools, 2nd ed.; Pearson/Addison-Wesley: Boston, 2007. (23) Rubin, N. C.; DePrince, A. E., III p† q: a tool for prototyping many-body methods for quantum chemistry. Mol. Phys. 2021, 119, e1954709. (24) Liebenthal, M. D.; Yuwono, S. H.; Koulias, L. N.; Li, R. R.; Rubin, N. C.; DePrince, A. E., III Automated Quantum Chemistry Code Generation with the p† q Package. J. Phys. Chem. A 2025, 129, 6679–6693. (25) Song, C.; Wang, L.-P.; Martı́nez, T. J. Automated Code Engine for Graphical Processing Units: Application to the Effective Core Potential Integrals and Gradients. J. Chem. Theory Comput. 2016, 12, 92–106. (26) Guerrero, R. D. RECURSUM: Automated Code Generation for Recurrence Relations Exceeds Expert Optimization via LayeredCodegen. arXiv:2604.18585 [cs.PL] 2026, (27) Hohenstein, E. G.; Parrish, R. M.; Martı́nez, T. J. Tensor hypercontraction density fitting. I. Quartic scaling second- and third-order Møller–Plesset perturbation theory. J. Chem. Phys. 2012, 137, 044103. (28) Feng, X.; Epifanovsky, E.; Gauss, J.; Krylov, A. I. Implementation of analytic gradients for CCSD and EOM-CCSD using Cholesky decomposition of the electron-repulsion
70
integrals and its application to disulfide radical anions. J. Chem. Phys. 2019, 151, 014110. (29) Almlöf, J.; Faegri, K.; Korsell, K. Principles for a direct SCF approach to LCAO–MO ab-initio calculations. J. Comput. Chem. 1982, 3, 385–399. (30) Warren, H. S., Jr. Hacker’s Delight, 2nd ed.; Addison-Wesley: Upper Saddle River, NJ, 2013. (31) Gouterman, M. Spectra of porphyrins. J. Mol. Spectrosc. 1961, 6, 138–163. (32) Rubio, M.; Roos, B. O.; Serrano-Andrés, L.; Merchán, M. Theoretical study of the electronic spectrum of magnesium-porphyrin. J. Chem. Phys. 1999, 110, 7202–7209. (33) Aikens, C. M.; Webb, S. P.; Bell, R. L.; Fletcher, G. D.; Schmidt, M. W.; Gordon, M. S. A derivation of the frozen-orbital unrestricted open-shell and restricted closed-shell second-order perturbation theory analytic gradient expressions. Theor. Chem. Acc. 2003, 110, 233–253.
71
TOC Graphic forward pass: densities
compile
d1
d2
F
∇Ek NACME
Goldstone term backward pass: relaxation ζ = transpose
A Goldstone term is a computation graph; its transpose—backpropagation—is the EOM-CC relaxation.
72