Notebooks · posts/casida-tdhf-and-tda-for-water

Excitation energies of water from the Casida equation, with and without de-excitations

Adam Zettel · September 25, 2026 · written with an agent

Water in the 6-31G basis, from its Hartree-Fock ground state to all 40 singlet and 40 triplet excitation energies, two ways. Time-dependent Hartree-Fock (the random-phase approximation) keeps the coupling $B$ between excitations and de-excitations; the Tamm-Dancoff approximation drops it and diagonalizes $A$ alone. The lowest singlet is at 0.344074 Ha (9.363 eV) in TDHF and 0.346161 Ha in TDA, and every one of the 160 energies matches PySCF's to $1 0^{- 8}$ Ha.

(A - B B - A) (X Y) = ω (X Y), (A - B)^{1/2} (A + B) (A - B)^{1/2} Z = ω^{2} Z

A molecule absorbs light at the energies of its excited states above the ground state, its excitation energies. Hartree-Fock gives only the ground state, a single set of occupied orbitals. Linear response asks how that state responds to a weak oscillating field: the response diverges at frequencies $ω$ that solve the eigenvalue problem above, the Casida equation. The unknowns are the amplitudes $X_{ia}$ for moving an electron from an occupied orbital $i$ to a virtual one $a$ , and $Y_{ia}$ for the reverse.

For a closed-shell molecule with real orbitals, orbital energies $ε$ , and two-electron integrals $(pq ∣ r s)$ in the molecular-orbital basis (chemists' notation),

A_{ia, j b} = δ_{ij} δ_{ab} (ε_{a} - ε_{i}) + 2 (ia ∣ j b) - (ij ∣ ab), B_{ia, j b} = 2 (ia ∣ j b) - (ib ∣ j a)

for singlets. For triplets the Coulomb terms $2 (ia ∣ j b)$ drop out, since a spin flip moves no net charge. With $B$ the equation is time-dependent Hartree-Fock (TDHF), also called the random-phase approximation (RPA). With $B = 0$ only $A$ is left, and its eigenvalues are the Tamm-Dancoff (TDA) energies, the same as configuration interaction with single excitations.

The program computes the integrals, the ground state and both spectra in this file, with lib/hf for the integrals and the self-consistent field. The reference is PySCF 2.5.0 (posts/refs/casida-tdhf-and-tda-for-water.py), which shares no code with it and solves both problems by Davidson iteration on matrix-vector products, without building $A$ or $B$ .

Look inside this program

Atlas
The program’s map: every function and which calls which.
Explore
Run it and step through what each part did, with the code, values and plots beside it.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 39–45

from hf.mole import molecule
from hf.molecule import mol_nelectron
from hf.scf import integrals
from hf.eri_store import es_unpack8
from hf.scf_loop import scf, scf_densities
from hf.opts import scf_defaults
from hf.stab import stab_fock_mo, stab_mo_eri, stab_orbitals

The ground state and its orbitals

The molecule is PySCF's water, written as PySCF takes it, and the 6-31G basis is found by its name in the table lib/basis/pople631g: 13 basis functions, 5 doubly occupied orbitals and 8 virtual ones, so $5 \times 8 = 40$ occupied-virtual pairs. The RHF energy matches PySCF's $- 75.98394849810573$ Ha. The orbitals are the eigenvectors of the converged Fock matrix, and the integrals $(pq ∣ r s)$ over them come from four quarter-transformations of the atomic-orbital integrals.

[2]lines 58–80

water = molecule("O 0 0 0; H 0 -0.757 0.587; H 0 0.757 0.587", "6-31g")
setup = integrals(water, scf_defaults())
H = setup.h
G = setup.eri
ground = scf(setup)
E_rhf = ground.energy
show E_rhf
nbf = len(H)
orb = stab_orbitals(setup, scf_densities(ground))
Fmo = stab_fock_mo(orb, H, G)
mo_eri = stab_mo_eri(es_unpack8(G, nbf), orb, 0, 0)
nocc = round(mol_nelectron(water) / 2)
nvir = nbf - nocc
nov = nocc * nvir
eps = zeros(nbf)
for p to nbf {
    eps[p] = Fmo[0, p, p]
}
show eps

test "the RHF energy is PySCF's -75.98394849810573 to 1e-9" {
    assert_near(E_rhf, -75.98394849810573, 1.0e-9)
}

Building $A$ and $B$ , and solving

The pair $ia$ is stored at index $i \cdot n_{vir} + (a - n_{occ})$ .

TDA is one symmetric eigenvalue problem, $A X = ω X$ .

For TDHF the $2 n \times 2 n$ problem is not symmetric, but it folds into a symmetric one of size $n$ . Adding and subtracting its two block rows gives $(A + B) (X + Y) = ω (X - Y)$ and $(A - B) (X - Y) = ω (X + Y)$ , so $(A - B) (A + B) (X + Y) = ω^{2} (X + Y)$ . When $A - B$ is positive definite, which holds when the Hartree-Fock state is stable against these rotations, it has a square root, and $Z = (A - B)^{- 1/2} (X + Y)$ turns the product into the symmetric matrix under the title. Its eigenvalues are $ω^{2}$ .

All 160 energies agree with PySCF's to $3.7 \times 1 0^{- 9}$ Ha, about what is left of the self-consistent field's convergence; the tests allow $1 0^{- 8}$ .

[3]lines 101–234

#| A and B for singlets (`spin` = 0) or triplets (`spin` = 1), written
#| into A and B from the orbital energies and the MO integrals.
def casida_ab(eps: Arr, M: Arr, nocc, nvir, spin, A: Arr, B: Arr) {
    coul = 2.0
    if spin > 0.5 {
        coul = 0.0
    }
    for i to nocc {
        for a to nvir {
            ia = i * nvir + a
            for j to nocc {
                for b to nvir {
                    jb = j * nvir + b
                    A[ia, jb] = coul * M[i, nocc + a, j, nocc + b] - M[i, j, nocc + a, nocc + b]
                    B[ia, jb] = coul * M[i, nocc + a, j, nocc + b] - M[i, nocc + b, j, nocc + a]
                }
            }
            A[ia, ia] += eps[nocc + a] - eps[i]
        }
    }
    return 0
}

#| C = X Y for square matrices.
def matmul_into(X: Arr, Y: Arr, C: Arr) {
    n = len(X)
    for i to n {
        for j to n {
            s = 0.0
            for k to n {
                s += X[i, k] * Y[k, j]
            }
            C[i, j] = s
        }
    }
    return 0
}

#| The TDA energies: the eigenvalues of A, ascending.
def tda_energies(A: Arr) {
    n = len(A)
    w = zeros(n)
    V = zeros(n, n)
    eig_sym(A + 0.0, w, V)
    sort(w)
    return w
}

#| The TDHF energies: omega = sqrt(lambda) for the eigenvalues lambda
#| of (A - B)^(1/2) (A + B) (A - B)^(1/2). Returns -1 in place of the
#| energies if A - B is not positive definite.
def tdhf_energies(A: Arr, B: Arr) {
    n = len(A)
    w = zeros(n)
    V = zeros(n, n)
    eig_sym(A - B, w, V)
    if min(w) <= 0.0 {
        return zeros(n) - 1.0
    }
    R = zeros(n, n)
    for i to n {
        for j to n {
            s = 0.0
            for k to n {
                s += V[i, k] * sqrt(w[k]) * V[j, k]
            }
            R[i, j] = s
        }
    }
    T = zeros(n, n)
    C = zeros(n, n)
    matmul_into(A + B, R, T)
    matmul_into(R, T, C)
    lam = zeros(n)
    eig_sym(C, lam, V)
    om = zeros(n)
    for k to n {
        om[k] = sqrt(lam[k])
    }
    sort(om)
    return om
}

A_s = zeros(nov, nov)
B_s = zeros(nov, nov)
casida_ab(eps, mo_eri, nocc, nvir, 0, A_s, B_s)
A_t = zeros(nov, nov)
B_t = zeros(nov, nov)
casida_ab(eps, mo_eri, nocc, nvir, 1, A_t, B_t)
s_tda = tda_energies(A_s)
s_tdhf = tdhf_energies(A_s, B_s)
t_tda = tda_energies(A_t)
t_tdhf = tdhf_energies(A_t, B_t)
HARTREE_EV = 27.211386245988
lowest_singlet_ev = s_tdhf[0] * HARTREE_EV
show nov
show lowest_singlet_ev

PYSCF_S_TDHF = [0.344073896628, 0.414613628694, 0.433066626728, 0.509323227106, 0.568680884883, 0.701941531252, 1.128806556941, 1.132160334289, 1.154582999530, 1.178856981854, 1.212201119944, 1.241207181902, 1.279014202889, 1.293065411055, 1.324682749963, 1.339695945976, 1.357645870338, 1.384802067573, 1.393741059918, 1.403238034836, 1.429962870047, 1.530829325635, 1.711259349922, 1.764341653599, 1.810672110269, 1.922282265948, 1.987891136260, 1.989412565866, 2.034512812717, 2.133877621188, 2.208729355547, 2.608180978545, 20.249465986572, 20.274171835043, 20.839432630354, 20.891124269144, 21.099768684428, 21.166188583007, 21.210294993184, 21.589655577140]
PYSCF_S_TDA = [0.346160646268, 0.417312694733, 0.436104114598, 0.512632037717, 0.570626213548, 0.711302432541, 1.131911718762, 1.132270713217, 1.157763254023, 1.182695980665, 1.213633893390, 1.242771446898, 1.279560761751, 1.294733958878, 1.326610798028, 1.341950087753, 1.361251575341, 1.393192516397, 1.396879533315, 1.404537102032, 1.447607649398, 1.537229483519, 1.711996305007, 1.770869880195, 1.873646212142, 1.923254131947, 1.989000360165, 1.990364251813, 2.042714277861, 2.156036557278, 2.209420590455, 2.643701258628, 20.249642780406, 20.274310681687, 20.839668376800, 20.891318787241, 21.099866977468, 21.166257883706, 21.210433865437, 21.589824806231]
PYSCF_T_TDHF = [0.306491899467, 0.366952973063, 0.389143731037, 0.430373420972, 0.504491208963, 0.555765151154, 0.954619810228, 1.072143598067, 1.098174184484, 1.118588097838, 1.140846393538, 1.145743271818, 1.149708340458, 1.200570903657, 1.210988703952, 1.232659852487, 1.283590409616, 1.308552012217, 1.325065797118, 1.338984895170, 1.377509897677, 1.418354815654, 1.510471355426, 1.601943812730, 1.714959134068, 1.853741267400, 1.933878845211, 1.943303504402, 1.954119606023, 2.071800470812, 2.164245337381, 2.480647938255, 20.184366056671, 20.210692019421, 20.745580010791, 20.808778741371, 21.052981725913, 21.152030016847, 21.160044471545, 21.542406352724]
PYSCF_T_TDA = [0.310920466443, 0.377683440260, 0.393748334010, 0.443103576263, 0.510118197001, 0.571637091366, 0.973670606684, 1.080367097903, 1.108669577126, 1.118673103703, 1.146403927520, 1.151514665188, 1.154708872451, 1.206975973268, 1.215838755385, 1.234665386598, 1.288448896581, 1.313879293377, 1.325791391161, 1.344333131897, 1.379797257831, 1.422319608884, 1.515325122561, 1.605987123814, 1.717494967357, 1.856144930100, 1.937842709741, 1.944020480719, 1.957038108998, 2.072922712202, 2.166120918081, 2.482182028729, 20.184459888354, 20.210792523016, 20.745665104182, 20.808865148557, 21.053033027246, 21.152041169783, 21.160113905331, 21.542459866451]

#| The largest |x_k - y_k|.
def max_gap(x: Arr, y: Arr) {
    g = 0.0
    for k to len(x) {
        g = max(g, abs(x[k] - y[k]))
    }
    return g
}

# the four lists together: every energy, and the worst disagreement with PySCF
n_energies = 4 * nov
worst_pyscf_gap = max(max(max_gap(s_tdhf, PYSCF_S_TDHF), max_gap(s_tda, PYSCF_S_TDA)), max(max_gap(t_tdhf, PYSCF_T_TDHF), max_gap(t_tda, PYSCF_T_TDA)))
show n_energies
show worst_pyscf_gap

test "all 160 energies agree with PySCF's to 3.7e-9 Ha" {
    assert_near(n_energies, 160, 0.0)
    assert_lt(worst_pyscf_gap, 3.7e-9)
}
test "all 40 singlet TDHF energies are PySCF's to 1e-8 Ha" {
    assert_lt(max_gap(s_tdhf, PYSCF_S_TDHF), 1.0e-8)
}
test "all 40 singlet TDA energies are PySCF's to 1e-8 Ha" {
    assert_lt(max_gap(s_tda, PYSCF_S_TDA), 1.0e-8)
}
test "all 40 triplet TDHF energies are PySCF's to 1e-8 Ha" {
    assert_lt(max_gap(t_tdhf, PYSCF_T_TDHF), 1.0e-8)
}
test "all 40 triplet TDA energies are PySCF's to 1e-8 Ha" {
    assert_lt(max_gap(t_tda, PYSCF_T_TDA), 1.0e-8)
}

#| A and B for singlets (`spin` = 0) or triplets (`spin` = 1), written
#| into A and B from the orbital energies and the MO integrals.
def casida_ab(eps: Arr, M: Arr, nocc, nvir, spin, A: Arr, B: Arr) {
    coul = 2.0
    if spin > 0.5 {
        coul = 0.0
    }
    for i to nocc {
        for a to nvir {
            ia = i * nvir + a
            for j to nocc {
                for b to nvir {
                    jb = j * nvir + b
                    A[ia, jb] = coul * M[i, nocc + a, j, nocc + b] - M[i, j, nocc + a, nocc + b]
                    B[ia, jb] = coul * M[i, nocc + a, j, nocc + b] - M[i, nocc + b, j, nocc + a]
                }
            }
            A[ia, ia] += eps[nocc + a] - eps[i]
        }
    }
    return 0
}

#| C = X Y for square matrices.
def matmul_into(X: Arr, Y: Arr, C: Arr) {
    n = len(X)
    for i to n {
        for j to n {
            s = 0.0
            for k to n {
                s += X[i, k] * Y[k, j]
            }
            C[i, j] = s
        }
    }
    return 0
}

#| The TDA energies: the eigenvalues of A, ascending.
def tda_energies(A: Arr) {
    n = len(A)
    w = zeros(n)
    V = zeros(n, n)
    eig_sym(A + 0.0, w, V)
    sort(w)
    return w
}

#| The TDHF energies: omega = sqrt(lambda) for the eigenvalues lambda
#| of (A - B)^(1/2) (A + B) (A - B)^(1/2). Returns -1 in place of the
#| energies if A - B is not positive definite.
def tdhf_energies(A: Arr, B: Arr) {
    n = len(A)
    w = zeros(n)
    V = zeros(n, n)
    eig_sym(A - B, w, V)
    if min(w) <= 0.0 {
        return zeros(n) - 1.0
    }
    R = zeros(n, n)
    for i to n {
        for j to n {
            s = 0.0
            for k to n {
                s += V[i, k] * sqrt(w[k]) * V[j, k]
            }
            R[i, j] = s
        }
    }
    T = zeros(n, n)
    C = zeros(n, n)
    matmul_into(A + B, R, T)
    matmul_into(R, T, C)
    lam = zeros(n)
    eig_sym(C, lam, V)
    om = zeros(n)
    for k to n {
        om[k] = sqrt(lam[k])
    }
    sort(om)
    return om
}

A_s = zeros(nov, nov)
B_s = zeros(nov, nov)
casida_ab(eps, mo_eri, nocc, nvir, 0, A_s, B_s)
A_t = zeros(nov, nov)
B_t = zeros(nov, nov)
casida_ab(eps, mo_eri, nocc, nvir, 1, A_t, B_t)
s_tda = tda_energies(A_s)
s_tdhf = tdhf_energies(A_s, B_s)
t_tda = tda_energies(A_t)
t_tdhf = tdhf_energies(A_t, B_t)
HARTREE_EV = 27.211386245988
lowest_singlet_ev = s_tdhf[0] * HARTREE_EV
show nov
show lowest_singlet_ev

PYSCF_S_TDHF = [0.344073896628, 0.414613628694, 0.433066626728, 0.509323227106, 0.568680884883, 0.701941531252, 1.128806556941, 1.132160334289, 1.154582999530, 1.178856981854, 1.212201119944, 1.241207181902, 1.279014202889, 1.293065411055, 1.324682749963, 1.339695945976, 1.357645870338, 1.384802067573, 1.393741059918, 1.403238034836, 1.429962870047, 1.530829325635, 1.711259349922, 1.764341653599, 1.810672110269, 1.922282265948, 1.987891136260, 1.989412565866, 2.034512812717, 2.133877621188, 2.208729355547, 2.608180978545, 20.249465986572, 20.274171835043, 20.839432630354, 20.891124269144, 21.099768684428, 21.166188583007, 21.210294993184, 21.589655577140]
PYSCF_S_TDA = [0.346160646268, 0.417312694733, 0.436104114598, 0.512632037717, 0.570626213548, 0.711302432541, 1.131911718762, 1.132270713217, 1.157763254023, 1.182695980665, 1.213633893390, 1.242771446898, 1.279560761751, 1.294733958878, 1.326610798028, 1.341950087753, 1.361251575341, 1.393192516397, 1.396879533315, 1.404537102032, 1.447607649398, 1.537229483519, 1.711996305007, 1.770869880195, 1.873646212142, 1.923254131947, 1.989000360165, 1.990364251813, 2.042714277861, 2.156036557278, 2.209420590455, 2.643701258628, 20.249642780406, 20.274310681687, 20.839668376800, 20.891318787241, 21.099866977468, 21.166257883706, 21.210433865437, 21.589824806231]
PYSCF_T_TDHF = [0.306491899467, 0.366952973063, 0.389143731037, 0.430373420972, 0.504491208963, 0.555765151154, 0.954619810228, 1.072143598067, 1.098174184484, 1.118588097838, 1.140846393538, 1.145743271818, 1.149708340458, 1.200570903657, 1.210988703952, 1.232659852487, 1.283590409616, 1.308552012217, 1.325065797118, 1.338984895170, 1.377509897677, 1.418354815654, 1.510471355426, 1.601943812730, 1.714959134068, 1.853741267400, 1.933878845211, 1.943303504402, 1.954119606023, 2.071800470812, 2.164245337381, 2.480647938255, 20.184366056671, 20.210692019421, 20.745580010791, 20.808778741371, 21.052981725913, 21.152030016847, 21.160044471545, 21.542406352724]
PYSCF_T_TDA = [0.310920466443, 0.377683440260, 0.393748334010, 0.443103576263, 0.510118197001, 0.571637091366, 0.973670606684, 1.080367097903, 1.108669577126, 1.118673103703, 1.146403927520, 1.151514665188, 1.154708872451, 1.206975973268, 1.215838755385, 1.234665386598, 1.288448896581, 1.313879293377, 1.325791391161, 1.344333131897, 1.379797257831, 1.422319608884, 1.515325122561, 1.605987123814, 1.717494967357, 1.856144930100, 1.937842709741, 1.944020480719, 1.957038108998, 2.072922712202, 2.166120918081, 2.482182028729, 20.184459888354, 20.210792523016, 20.745665104182, 20.808865148557, 21.053033027246, 21.152041169783, 21.160113905331, 21.542459866451]

#| The largest |x_k - y_k|.
def max_gap(x: Arr, y: Arr) {
    g = 0.0
    for k to len(x) {
        g = max(g, abs(x[k] - y[k]))
    }
    return g
}

# the four lists together: every energy, and the worst disagreement with PySCF
n_energies = 4 * nov
worst_pyscf_gap = max(max(max_gap(s_tdhf, PYSCF_S_TDHF), max_gap(s_tda, PYSCF_S_TDA)), max(max_gap(t_tdhf, PYSCF_T_TDHF), max_gap(t_tda, PYSCF_T_TDA)))
show n_energies
show worst_pyscf_gap

test "all 160 energies agree with PySCF's to 3.7e-9 Ha" {
    assert_near(n_energies, 160, 0.0)
    assert_lt(worst_pyscf_gap, 3.7e-9)
}
test "all 40 singlet TDHF energies are PySCF's to 1e-8 Ha" {
    assert_lt(max_gap(s_tdhf, PYSCF_S_TDHF), 1.0e-8)
}
test "all 40 singlet TDA energies are PySCF's to 1e-8 Ha" {
    assert_lt(max_gap(s_tda, PYSCF_S_TDA), 1.0e-8)
}
test "all 40 triplet TDHF energies are PySCF's to 1e-8 Ha" {
    assert_lt(max_gap(t_tdhf, PYSCF_T_TDHF), 1.0e-8)
}
test "all 40 triplet TDA energies are PySCF's to 1e-8 Ha" {
    assert_lt(max_gap(t_tda, PYSCF_T_TDA), 1.0e-8)
}

What $B$ does

$B$ couples each excitation to the de-excitations, and its effect is to pull the energies down. For the lowest root this always holds. With $Λ = (A B B A)$ , the lowest TDHF $ω$ is the minimum of $(X, Y)^{T} Λ (X, Y) / (X^{T} X - Y^{T} Y)$ over vectors with $X^{T} X > Y^{T} Y$ ; the vectors with $Y = 0$ are among them, and over those the ratio is $X^{T} A X / X^{T} X$ , whose minimum is the lowest TDA energy. For water every root comes out lower in TDHF than in TDA, by 0.07 to 63 millihartree for the singlets and 0.01 to 19 for the triplets (posts/refs prints the same ranges from PySCF); the largest is the 25th singlet's, 63 millihartree. The 8 roots near 20 Ha, which move an electron out of the oxygen 1s orbital, shift by less than 0.25 millihartree.

The triplets sit below the singlets. The singlet $A$ and $B$ carry the extra term $2 (ia ∣ j b)$ , whose diagonal $2 (ia ∣ ia)$ is the Coulomb self-repulsion of the charge distribution $ϕ_{i} ϕ_{a}$ and is positive. The lowest triplet is 0.3065 Ha, 1.023 eV below the lowest singlet.

The first figure shows the 32 valence excitation energies of each kind. The second shows, root by root, how far TDA lies above TDHF.

[4]lines 259–357

s_shift = s_tda - s_tdhf
# the roots below the 8 oxygen 1s excitations near 20 Ha are the valence ones
N_VAL = nov - nvir
t_shift = t_tda - t_tdhf
# the same shifts in millihartree, as the prose states them, and the largest
# shift of the 8 oxygen 1s roots of either kind
shift_singlet_min_mh = 1000.0 * min(s_shift)
shift_singlet_max_mh = 1000.0 * max(s_shift)
shift_triplet_min_mh = 1000.0 * min(t_shift)
shift_triplet_max_mh = 1000.0 * max(t_shift)
core_shift_max_mh = 0.0
for k from N_VAL to nov {
    core_shift_max_mh = max(core_shift_max_mh, 1000.0 * max(s_shift[k], t_shift[k]))
}
largest_shift_root = 0
for k to nov {
    if s_shift[k] > s_shift[largest_shift_root] {
        largest_shift_root = k
    }
}
show shift_singlet_min_mh
show shift_singlet_max_mh
show shift_triplet_min_mh
show shift_triplet_max_mh
show core_shift_max_mh
show N_VAL
shift_singlet_min = min(s_shift)
shift_singlet_max = max(s_shift)
shift_triplet_min = min(t_shift)
shift_triplet_max = max(t_shift)
show shift_singlet_min
show shift_singlet_max
show shift_triplet_min
show shift_triplet_max
singlet_triplet_gap_ev = (s_tdhf[0] - t_tdhf[0]) * HARTREE_EV
show singlet_triplet_gap_ev

test "every TDHF root is below its TDA root: by 6.930e-5 to 6.297e-2 Ha for singlets, 1.115e-5 to 1.905e-2 for triplets" {
    assert_near(min(s_shift), 6.930e-5, 1.0e-7)
    assert_near(max(s_shift), 6.297e-2, 1.0e-5)
    assert_near(min(t_shift), 1.115e-5, 1.0e-7)
    assert_near(max(t_shift), 1.905e-2, 1.0e-5)
}
test "TDA lies above TDHF by 0.07 to 63 millihartree for singlets, largest at the 25th, and 0.01 to 19 for triplets" {
    assert_near(shift_singlet_min_mh, 0.0693, 1.0e-4)
    assert_near(shift_singlet_max_mh, 62.97, 0.01)
    assert_near(largest_shift_root + 1, 25, 0.0)
    assert_near(shift_triplet_min_mh, 0.0112, 1.0e-4)
    assert_near(shift_triplet_max_mh, 19.05, 0.01)
}
test "32 of the 40 roots are valence, below the 8 oxygen 1s roots near 20 Ha" {
    assert_near(N_VAL, 32, 0.0)
    assert_lt(s_tdhf[N_VAL - 1], 3.0)
    assert_gt(s_tdhf[N_VAL], 20.0)
}
test "the 8 oxygen 1s excitations move by less than 0.25 millihartree" {
    assert_lt(core_shift_max_mh, 0.25)
    for k from N_VAL to nov {
        assert_lt(s_shift[k], 2.5e-4)
        assert_lt(t_shift[k], 2.5e-4)
    }
}
test "the lowest singlet is 9.3627 eV in TDHF, and the lowest triplet 1.0227 eV below it" {
    assert_near(lowest_singlet_ev, 0.344073896628 * HARTREE_EV, 1.0e-6)
    assert_near(singlet_triplet_gap_ev, (0.344073896628 - 0.306491899467) * HARTREE_EV, 1.0e-6)
}

# the 32 valence roots (below the oxygen 1s excitations near 20 Ha), in eV
idx = zeros(N_VAL)
s_rpa_ev = zeros(N_VAL)
s_tda_ev = zeros(N_VAL)
t_rpa_ev = zeros(N_VAL)
t_tda_ev = zeros(N_VAL)
s_d = zeros(N_VAL)
t_d = zeros(N_VAL)
for k to N_VAL {
    idx[k] = k + 1
    s_rpa_ev[k] = s_tdhf[k] * HARTREE_EV
    s_tda_ev[k] = s_tda[k] * HARTREE_EV
    t_rpa_ev[k] = t_tdhf[k] * HARTREE_EV
    t_tda_ev[k] = t_tda[k] * HARTREE_EV
    s_d[k] = 1000.0 * s_shift[k]
    t_d[k] = 1000.0 * t_shift[k]
}
figure("Water, 6-31G: valence excitation energies", "root, in order of energy", "excitation energy (eV)")
plot(idx, s_rpa_ev, "singlet, TDHF")
plot(idx, s_tda_ev, "singlet, TDA")
plot(idx, t_rpa_ev, "triplet, TDHF")
plot(idx, t_tda_ev, "triplet, TDA")
style("singlet, TDA", "dash")
style("triplet, TDA", "dash")

figure("How far TDA lies above TDHF, root by root", "root, in order of energy", "TDA - TDHF (millihartree)")
yscale("log")
plot(idx, s_d, "singlets")
plot(idx, t_d, "triplets")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

s_shift = s_tda - s_tdhf
# the roots below the 8 oxygen 1s excitations near 20 Ha are the valence ones
N_VAL = nov - nvir
t_shift = t_tda - t_tdhf
# the same shifts in millihartree, as the prose states them, and the largest
# shift of the 8 oxygen 1s roots of either kind
shift_singlet_min_mh = 1000.0 * min(s_shift)
shift_singlet_max_mh = 1000.0 * max(s_shift)
shift_triplet_min_mh = 1000.0 * min(t_shift)
shift_triplet_max_mh = 1000.0 * max(t_shift)
core_shift_max_mh = 0.0
for k from N_VAL to nov {
    core_shift_max_mh = max(core_shift_max_mh, 1000.0 * max(s_shift[k], t_shift[k]))
}
largest_shift_root = 0
for k to nov {
    if s_shift[k] > s_shift[largest_shift_root] {
        largest_shift_root = k
    }
}
show shift_singlet_min_mh
show shift_singlet_max_mh
show shift_triplet_min_mh
show shift_triplet_max_mh
show core_shift_max_mh
show N_VAL
shift_singlet_min = min(s_shift)
shift_singlet_max = max(s_shift)
shift_triplet_min = min(t_shift)
shift_triplet_max = max(t_shift)
show shift_singlet_min
show shift_singlet_max
show shift_triplet_min
show shift_triplet_max
singlet_triplet_gap_ev = (s_tdhf[0] - t_tdhf[0]) * HARTREE_EV
show singlet_triplet_gap_ev

test "every TDHF root is below its TDA root: by 6.930e-5 to 6.297e-2 Ha for singlets, 1.115e-5 to 1.905e-2 for triplets" {
    assert_near(min(s_shift), 6.930e-5, 1.0e-7)
    assert_near(max(s_shift), 6.297e-2, 1.0e-5)
    assert_near(min(t_shift), 1.115e-5, 1.0e-7)
    assert_near(max(t_shift), 1.905e-2, 1.0e-5)
}
test "TDA lies above TDHF by 0.07 to 63 millihartree for singlets, largest at the 25th, and 0.01 to 19 for triplets" {
    assert_near(shift_singlet_min_mh, 0.0693, 1.0e-4)
    assert_near(shift_singlet_max_mh, 62.97, 0.01)
    assert_near(largest_shift_root + 1, 25, 0.0)
    assert_near(shift_triplet_min_mh, 0.0112, 1.0e-4)
    assert_near(shift_triplet_max_mh, 19.05, 0.01)
}
test "32 of the 40 roots are valence, below the 8 oxygen 1s roots near 20 Ha" {
    assert_near(N_VAL, 32, 0.0)
    assert_lt(s_tdhf[N_VAL - 1], 3.0)
    assert_gt(s_tdhf[N_VAL], 20.0)
}
test "the 8 oxygen 1s excitations move by less than 0.25 millihartree" {
    assert_lt(core_shift_max_mh, 0.25)
    for k from N_VAL to nov {
        assert_lt(s_shift[k], 2.5e-4)
        assert_lt(t_shift[k], 2.5e-4)
    }
}
test "the lowest singlet is 9.3627 eV in TDHF, and the lowest triplet 1.0227 eV below it" {
    assert_near(lowest_singlet_ev, 0.344073896628 * HARTREE_EV, 1.0e-6)
    assert_near(singlet_triplet_gap_ev, (0.344073896628 - 0.306491899467) * HARTREE_EV, 1.0e-6)
}

# the 32 valence roots (below the oxygen 1s excitations near 20 Ha), in eV
idx = zeros(N_VAL)
s_rpa_ev = zeros(N_VAL)
s_tda_ev = zeros(N_VAL)
t_rpa_ev = zeros(N_VAL)
t_tda_ev = zeros(N_VAL)
s_d = zeros(N_VAL)
t_d = zeros(N_VAL)
for k to N_VAL {
    idx[k] = k + 1
    s_rpa_ev[k] = s_tdhf[k] * HARTREE_EV
    s_tda_ev[k] = s_tda[k] * HARTREE_EV
    t_rpa_ev[k] = t_tdhf[k] * HARTREE_EV
    t_tda_ev[k] = t_tda[k] * HARTREE_EV
    s_d[k] = 1000.0 * s_shift[k]
    t_d[k] = 1000.0 * t_shift[k]
}
figure("Water, 6-31G: valence excitation energies", "root, in order of energy", "excitation energy (eV)")
plot(idx, s_rpa_ev, "singlet, TDHF")
plot(idx, s_tda_ev, "singlet, TDA")
plot(idx, t_rpa_ev, "triplet, TDHF")
plot(idx, t_tda_ev, "triplet, TDA")
style("singlet, TDA", "dash")
style("triplet, TDA", "dash")

figure("How far TDA lies above TDHF, root by root", "root, in order of energy", "TDA - TDHF (millihartree)")
yscale("log")
plot(idx, s_d, "singlets")
plot(idx, t_d, "triplets")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

The ground state and its orbitals

Building A and B, and solving

What B does

Building $A$ and $B$ , and solving

What $B$ does