feat(cis): Davidson iterative eigensolver wired into CIS / TDA

Tier 3 follow-up. CIS gains `useDavidson: true` opt-in that avoids
building the full A matrix on the OV space. Closes the
"Davidson / Krylov eigensolver" row in LIMITATIONS.md — Davidson
is now available across **every** eigenvalue-method path:
EE/IP/EA-EOM-CCSD (commits 6b9a96f, a35baeb) + CIS / TDA (this).

Implementation:
- New `davidsonExtract` path in `runCIS`. Builds a matvec from the
  same `eri()` accessor used by `buildA`:
    (Av)_{ia} = (ε_a − ε_i) v_ia
              + Σ_{jb} [s_Coul · (ai|bj) − (ij|ab)] v_jb
  s_Coul = 2 for singlet, 0 for triplet (same as dense path).
- Diagonal preconditioner is the Koopmans gap (ε_a − ε_i).
- Default off (`useDavidson` undefined) → existing dense path,
  bit-identical to before.

Verified by `cis-davidson.test.ts` on H₂O HF/STO-3G:
- k=3 singlets: Davidson vs dense agree to < 1e-6 Ha
- k=3 triplets: same

Also a co-located type-fix in `cphf.ts` (caught by the typecheck):
the symmetrize loop's array indexing was triggering
`noUncheckedIndexedAccess`; cleaned up with hoisted refs.

Full suite: 446 passing / 1 skipped (+2 from this commit), no
regressions.

Tier 3 progress: 13 items shipped this session.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: abgnydn

LIMITATIONS.md

@@ -124,7 +124,7 @@ will silently truncate large dispatches.
 | ccData / QC-Schema compatibility | Output not consumable by external tooling | Tier 3 |
 | FAIR / Zenodo DOIs per release | Citations point at GitHub tag, not DOI | Tier 3 — set up CI workflow |
 | Aux-basis density fitting | We have CD-DF but not JKFIT / RIFIT integral path | Tier 3 — needs 3-index ERI routine |
-| Davidson / Krylov eigensolver | OK for n_occ·n_virt ≤ ~500; dense above. Measured 2026-05: EOM-CCSD on CH₂O / HCN at STO-3G (dim 3488 / 2660) takes 15-30 min per molecule on M2 Pro because the Hessenberg + Wilkinson QR is dense. Davidson would make these ~minutes. | Tier 3 · high impact |
+| Davidson / Krylov eigensolver | **Shipped 2026-05**. Block Davidson (`src/manybody/davidson.ts`) wired into EE/IP/EA-EOM-CCSD (commits 6b9a96f, a35baeb) and now also CIS / TDA (this commit) behind `useDavidson: true`. Dense + Davidson agree on k roots to ≤ 1e-6 Ha across all paths. Unblocks large-basis CIS/TDDFT + bigger EOM-CCSD systems. | shipped ✓ |
 | Multi-node parallel | One tab, one GPU | Tier 4 — substrate is Phase D WebRTC |
 | Periodic boundary conditions | No solids, no surfaces | Tier 4 |
 | Spin-orbit coupling / X2C / DKH | No heavy elements | Tier 4 |

src/chemistry/cis.ts

@@ -32,6 +32,7 @@
 import type { MolecularIntegrals } from "./cg-molecular.js";
 import { transformERIToMO } from "./mp2.js";
 import { eigsymmetric } from "../manybody/dense-eig.js";
+import { davidson } from "../manybody/davidson.js";
 
 export interface HFLike {
   /** AO → MO coefficients (row-major, column i = MO i). */
@@ -50,6 +51,13 @@ export interface CISOpts {
   readonly nRoots?: number;
   /** Spin sector (default both). Singlet only is the common case. */
   readonly spin?: "singlet" | "triplet" | "both";
+  /** Use the block Davidson iterative eigensolver instead of building
+   *  the full A matrix and calling `eigsymmetric`. Recommended for
+   *  dim ≳ 200 (cc-pVDZ on multi-atom systems). Default false
+   *  preserves the existing dense behavior. */
+  readonly useDavidson?: boolean;
+  /** Davidson residual-norm tolerance. Default 1e-7. */
+  readonly davidsonTol?: number;
 }
 
 export interface CISStateBlock {
@@ -146,9 +154,55 @@ export function runCIS(
     return { energies, amplitudes };
   };
 
+  /** Davidson path — never materializes A. Builds a matvec from
+   *  the CIS A-block formula and supplies the Koopmans diagonal
+   *  preconditioner (ε_a − ε_i). */
+  const davidsonExtract = (sCoul: number): CISStateBlock => {
+    const eps = hf.orbitalEnergies;
+    const diagonal = new Float64Array(dim);
+    for (let i = 0; i < nOcc; i++) {
+      for (let a = 0; a < nVirt; a++) {
+        const ai = nOcc + a;
+        diagonal[i * nVirt + a] = eps[ai]! - eps[i]!;
+      }
+    }
+    const matvec = (v: Float64Array): Float64Array => {
+      // (Av)_{ia} = (ε_a − ε_i) v_ia
+      //           + Σ_{jb} [sCoul · (ia|jb) − (ij|ab)] v_jb
+      const out = new Float64Array(dim);
+      for (let i = 0; i < nOcc; i++) {
+        for (let a = 0; a < nVirt; a++) {
+          const ai = nOcc + a;
+          const eOrb = eps[ai]! - eps[i]!;
+          let s = eOrb * v[i * nVirt + a]!;
+          for (let j = 0; j < nOcc; j++) {
+            for (let b = 0; b < nVirt; b++) {
+              const bj = nOcc + b;
+              const coul = sCoul !== 0 ? sCoul * eri(i, ai, j, bj) : 0;
+              const exch = eri(i, j, ai, bj);
+              s += (coul - exch) * v[j * nVirt + b]!;
+            }
+          }
+          out[i * nVirt + a] = s;
+        }
+      }
+      return out;
+    };
+    const dav = davidson(dim, matvec, diagonal, {
+      k: nRoots,
+      tol: opts.davidsonTol ?? 1e-7,
+      maxIter: 200,
+      filter: (_re, im) => Math.abs(im) < 1e-6,
+    });
+    return { energies: dav.energies, amplitudes: dav.vectors };
+  };
+
   const empty: CISStateBlock = { energies: new Float64Array(0), amplitudes: new Float64Array(0) };
-  const singlet = (spin === "singlet" || spin === "both") ? diagAndExtract(buildA(2)) : empty;
-  const triplet = (spin === "triplet" || spin === "both") ? diagAndExtract(buildA(0)) : empty;
+  const extract = opts.useDavidson
+    ? davidsonExtract
+    : (sCoul: number): CISStateBlock => diagAndExtract(buildA(sCoul));
+  const singlet = (spin === "singlet" || spin === "both") ? extract(2) : empty;
+  const triplet = (spin === "triplet" || spin === "both") ? extract(0) : empty;
 
   return { nOccupied: nOcc, nVirtual: nVirt, singlet, triplet };
 }

src/chemistry/cphf.ts

@@ -187,7 +187,9 @@ export function cphfPolarizability(
   for (let x = 0; x < 3; x++) {
     for (let y = 0; y < 3; y++) {
       let s = 0;
-      for (let k = 0; k < dim; k++) s += X[x][k]! * muOV[y][k]!;
+      const Xx = X[x]!;
+      const Fy = muOV[y]!;
+      for (let k = 0; k < dim; k++) s += Xx[k]! * Fy[k]!;
       alpha[x * 3 + y] = 4 * s;
     }
   }

tests/chemistry/cis-davidson.test.ts

@@ -0,0 +1,61 @@
+// CIS Davidson eigensolver cross-check.
+//
+// Verifies that the iterative `useDavidson: true` path returns the
+// same lowest singlet + triplet excitation energies as the dense
+// `eigsymmetric` path. Davidson avoids building the full A matrix,
+// which scales as O(dim²) memory — important for cc-pVDZ on
+// multi-atom systems where dim grows past a few hundred.
+
+import { describe, expect, test } from "vitest";
+import { computeMolecularIntegrals } from "../../src/chemistry/cg-molecular.js";
+import { moleculeToShellsNuclei, type Atom } from "../../src/chemistry/atoms.js";
+import { runRHFSCF } from "../../src/chemistry/hf-scf.js";
+import { runCIS } from "../../src/chemistry/cis.js";
+
+describe("CIS — Davidson vs dense", () => {
+  test("H₂O HF/STO-3G — k=3 singlets agree to ≤ 1e-7 Ha", () => {
+    const half = (104.52 / 2) * Math.PI / 180;
+    const xH = 0.9572 * Math.sin(half);
+    const zH = 0.9572 * Math.cos(half);
+    const atoms: Atom[] = [
+      { symbol: "O", pos: [0, 0, 0] },
+      { symbol: "H", pos: [ xH, 0, zH] },
+      { symbol: "H", pos: [-xH, 0, zH] },
+    ];
+    const { shells, nuclei, nElectrons } = moleculeToShellsNuclei(atoms);
+    const integrals = computeMolecularIntegrals(shells, nuclei);
+    const hf = runRHFSCF(integrals, nElectrons, {
+      useDIIS: true, maxIter: 200, energyTol: 1e-10, densityTol: 1e-8,
+    });
+
+    const dense = runCIS(integrals, hf, { nRoots: 3, spin: "singlet" });
+    const dav   = runCIS(integrals, hf, { nRoots: 3, spin: "singlet", useDavidson: true, davidsonTol: 1e-9 });
+
+    expect(dense.singlet.energies.length).toBe(3);
+    expect(dav.singlet.energies.length).toBe(3);
+    for (let k = 0; k < 3; k++) {
+      expect(Math.abs(dav.singlet.energies[k]! - dense.singlet.energies[k]!)).toBeLessThan(1e-6);
+    }
+  }, 60_000);
+
+  test("H₂O HF/STO-3G — k=3 triplets agree to ≤ 1e-6 Ha", () => {
+    const half = (104.52 / 2) * Math.PI / 180;
+    const xH = 0.9572 * Math.sin(half);
+    const zH = 0.9572 * Math.cos(half);
+    const atoms: Atom[] = [
+      { symbol: "O", pos: [0, 0, 0] },
+      { symbol: "H", pos: [ xH, 0, zH] },
+      { symbol: "H", pos: [-xH, 0, zH] },
+    ];
+    const { shells, nuclei, nElectrons } = moleculeToShellsNuclei(atoms);
+    const integrals = computeMolecularIntegrals(shells, nuclei);
+    const hf = runRHFSCF(integrals, nElectrons, {
+      useDIIS: true, maxIter: 200, energyTol: 1e-10, densityTol: 1e-8,
+    });
+    const dense = runCIS(integrals, hf, { nRoots: 3, spin: "triplet" });
+    const dav   = runCIS(integrals, hf, { nRoots: 3, spin: "triplet", useDavidson: true });
+    for (let k = 0; k < 3; k++) {
+      expect(Math.abs(dav.triplet.energies[k]! - dense.triplet.energies[k]!)).toBeLessThan(1e-6);
+    }
+  }, 60_000);
+});