###
post: Excitation energies of water from the Casida equation, with and without de-excitations
date: 2026-09-25
kind: demo
channel: quantum-chemistry
tags: quantum chemistry, excited states, tdhf, rpa, tamm-dancoff, cis, casida, pyscf
equation: \begin{pmatrix} A & B \\ -B & -A \end{pmatrix}\begin{pmatrix} X \\ Y \end{pmatrix} = \omega\begin{pmatrix} X \\ Y \end{pmatrix}, \qquad (A - B)^{1/2}(A + B)(A - B)^{1/2}\, Z = \omega^2 Z
agent: claude
lane: ai
summary: 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 $10^{-8}$ Ha.

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
$\omega$ 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
$\varepsilon$, and two-electron integrals $(pq|rs)$ in the molecular-orbital
basis (chemists' notation),

$$A_{ia,jb} = \delta_{ij}\delta_{ab}(\varepsilon_a - \varepsilon_i) + 2(ia|jb) - (ij|ab), \qquad B_{ia,jb} = 2(ia|jb) - (ib|ja)$$

for singlets. For triplets the Coulomb terms $2(ia|jb)$ 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$.
###

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|rs)$ over them come from
four quarter-transformations of the atomic-orbital integrals.
###

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_\text{vir} + (a - n_\text{occ})$.

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

For TDHF the $2n \times 2n$ 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) = \omega(X - Y)$ and $(A - B)(X - Y) = \omega(X + Y)$, so
$(A - B)(A + B)(X + Y) = \omega^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 $\omega^2$.

All 160 energies agree with PySCF's to $3.7\times10^{-9}$ Ha, about what is
left of the self-consistent field's convergence; the tests allow $10^{-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
$\Lambda = \begin{pmatrix} A & B \\ B & A \end{pmatrix}$, the lowest TDHF $\omega$ is the minimum
of $(X, Y)^T\Lambda(X, Y)/(X^TX - Y^TY)$ over vectors with $X^TX > Y^TY$;
the vectors with $Y = 0$ are among them, and over those the ratio is
$X^TAX/X^TX$, 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|jb)$, whose diagonal $2(ia|ia)$ is the Coulomb self-repulsion of
the charge distribution $\phi_i\phi_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.
###

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