Casida's Equation

Casida's equation gives you a way to calculate the excitation energies and oscillator strengths from a DFT calculation. You know how when you push a swing, you have to push at the right frequency or the amplitude doesn't build? Instead of pushing a swing or a pendulum, Casida's equation gives you a way to push the electron cloud of the system. The twist is that the cloud is a quantum system — there are only certain photon frequencies at which it responds at all, and when one of those hits, the response isn't a gradual buildup like a swing; the cloud jumps discretely to an excited state. Those photon frequencies are the excitation energies.¹ How strongly the field couples to each one is the oscillator strength. Casida's equation is the matrix problem that takes the ground-state Kohn-Sham orbitals as input and computes both.

Setup: linear response of a closed-shell reference

Setup: a molecule in its ground state, solved by the standard SCF machinery (Kohn-Sham for DFT, Hartree-Fock for HF). You get a stack of orbitals. The bottom $N/2$ are occupied by paired electrons (${\psi }_i$, energy ${\epsilon }_i$); the rest are empty (virtuals ${\psi }_a$, energy ${\epsilon }_a$). That's the static ground state — what a normal DFT calculation gives you.

Quick aside on virtual orbitals: they're real orbitals — don't let the occupied orbitals hear that, but yeah, they are. When you solve an SCF problem you get a set of orbitals; some get filled with electrons, some don't, and the empty ones are called "virtual." Think of them as extra seats in a room: they exist, they're available, no one's sitting there right now. Why nobody's in them: promoting an electron from occupied $i$ into virtual $a$ costs energy ${\epsilon }_a-{\epsilon }_i$, and in the ground state there's nothing to pay that bill — no field, no thermal kick, nothing. So the electrons fill the bottom of the stack and the rest stay empty. Exciting the molecule means supplying exactly that bill from outside — which is what the oscillating field is going to do in the next paragraph.

Now turn on a weak oscillating field — a laser shining on the molecule, say. The electrons respond. To first order in the field strength, the response at each driving frequency $\omega$ splits into two pieces. $X_{ia}$ tells you how strongly the field is mixing occupied orbital $i$ with virtual orbital $a$ — that's an electron getting promoted. $Y_{ia}$ tells you the time-reverse: how strongly the same $i\leftrightarrow a$ coupling shows up in the opposite direction, an electron coming back down. Both appear together at every frequency; linear-response quantum mechanics is symmetric under time reversal, so you can't have one without the other.

Plug this into the time-dependent Kohn-Sham equation and keep only first-order terms in the field. The algebra — worked out in the deep dive at the bottom of the page — gives a coupled matrix equation for $X$ and $Y$. In block form:

What in the name of Charles Hermite is going on here? Why is this matrix non-Hermitian? Look closely. The matrix on the left is actually Hermitian. The culprit is one entry, hiding in plain sight: the -1 in the lower-right of the metric on the right-hand side. Without that minus sign the equation would be an ordinary Hermitian eigenvalue problem. With it, the metric is indefinite (eigenvalues $+1$ and $-1$ rather than both $+1$), and folding it into the left-hand matrix to turn this into an ordinary eigenvalue problem produces something non-Hermitian. The minus sign isn't an accident; it encodes time-reversal symmetry. Every positive-frequency root $\omega$ comes with a $-\omega$ partner representing the same transition run backwards. For real orbitals (the usual case) $A$ and $B$ are real symmetric, and:

Two pieces to unpack. The diagonal of $A$ is the orbital-energy gap — what it would cost to promote an electron from $i$ to $a$ if you ignored interactions. Everything else, the $(ia\Vert jb)$ terms, is the interaction part: how the $i\to a$ transition density Coulomb-repels and exchanges with the $j\to b$ transition density. TDHF uses Coulomb minus exchange; TDDFT uses Coulomb plus the XC kernel $f_{xc}$. The two blocks differ in what they couple. $A$ couples excitations to other excitations. $B$ couples excitations to de-excitations — the time-reversal partner of the same physics — and the Tamm-Dancoff approximation just throws that block away.

Spin adaptation: singlets and triplets

A closed-shell molecule has equal numbers of $\alpha$ and $\beta$ electrons paired in every orbital — full spin symmetry. That symmetry makes the Casida matrix split into independent blocks. The blocks correspond to singlet excitations (spin-symmetric: same amplitude for $\alpha \to \alpha$ and $\beta \to \beta$ transitions) and triplet excitations (antisymmetric: opposite-sign amplitudes for $\alpha \to \alpha$ and $\beta \to \beta$, plus two more components that flip spin). The three triplets are degenerate by symmetry, so you solve them all by solving one. Net: two independent problems of size $N_{\text{occ}}{N}_{\text{vir}}$ each, instead of one giant problem four times the size:

The factor of 2 in front of $(ia|jb)$ on the singlet side comes from the spin-symmetric combination; the triplet has no Coulomb piece in $A$ at all (the symmetric Coulomb interaction cancels between the $\alpha \to \alpha$ and $\beta \to \beta$ components of a triplet). What's left of the triplet is purely exchange — this is the origin of Hund's rule in the orbital language: triplets sit lower than singlets because they lose Coulomb repulsion.

Why it's a pseudo-eigenvalue problem

The minus sign in the metric is awkward, but when $A-B$ is positive-definite — which it is whenever the ground state is a stable minimum and not a saddle — you can fold the equation into an ordinary Hermitian eigenvalue problem half the size. Here is the algebra.

Write the block equation out as two coupled vector equations. With $A$ and $B$ real symmetric:

Add the two equations and subtract them. The sum gives one equation in $X+Y$ and $X-Y$, the difference gives the other:

Solve the second for $X+Y$ (using that $A-B$ is invertible) and substitute into the first. The cross-coupling drops out and what's left is an ordinary eigenvalue problem for $X-Y$ alone, with eigenvalues ${\omega }^2$:

The operator $(A+B)(A-B)$ is not Hermitian — products of Hermitian operators usually aren't. But it is similar to a Hermitian operator. Take the symmetric square root $(A-B{)}^{1/2}$ (which exists and is real because $A-B$ is positive-definite), define a rescaled eigenvector

and conjugate the eigenvalue problem by $(A-B{)}^{1/2}$. The result is Hermitian:

The eigenvalues are ${\omega }^2$, so excitation energies come out as square roots. The eigenvector $Z$ recovers the physical amplitudes via $X-Y=(A-B{)}^{-1/2}Z$ and $X+Y={\omega }^{-1}(A-B{)}^{1/2}Z$. A triplet instability — an eigenvalue of $A-B$ going negative — is the signature that the closed-shell reference is unstable against a spin-symmetry-breaking deformation, and Casida loudly tells you so by failing to take a real square root.

The Tamm-Dancoff approximation

Drop $B$ and the whole problem collapses to ordinary Hermitian diagonalization of $A$:

This is the Tamm-Dancoff approximation (TDA), or "configuration-interaction singles with the DFT orbital basis" depending on who's talking. Cheaper, more numerically stable (no square-root of a possibly-near-singular matrix), and immune to the triplet-instability failure mode. The price: TDA excitation energies are systematically too high, because the de-excitation coupling $B$ would have lowered the levels through a second-order-like correction and that correction is gone. For valence singlets the TDA error is a few tenths of an eV; for charge-transfer states and Rydberg-like states it can be larger.

Oscillator strengths and transition densities

Once you've got a Casida root $({\omega }_n,X_n,Y_n)$, the transition density — the spatial pattern of where the electron sloshes during this excitation — is a weighted sum over occupied-virtual orbital products:

Multiply this by $\mathbf{r}$ and integrate to get the transition dipole — that's what couples to photons. The oscillator strength $f_n=\frac{2}{3}{\omega }_n|⟨0|\mathbf{r}|n⟩{|}^2$ packages it into a dimensionless intensity. The absorption spectrum is then one delta function per root, weighted by $f_n$ and parked at ${\omega }_n$; broaden them into Lorentzians and you've got something to lay over a measured spectrum.

H₂ in a minimal basis: closed form

Before turning on PySCF, it's worth doing the smallest possible case by hand. H₂ in a minimal basis — one s-function per atom, two basis functions total — collapses the entire Casida equation to one number per spin block. The result is closed form, and it makes explicit where the singlet-triplet splitting comes from. This is the case Casida worked through in the original 1995 paper.

Symmetry forces the molecular orbitals to be the bonding and antibonding combinations of the two s-functions:

One occupied orbital (${\psi }_g$, bonding), one virtual orbital (${\psi }_u$, antibonding). The Casida problem dimension per spin sector is $N_{\mathrm{occ}}\cdot N_{\mathrm{vir}}=1$ — every matrix in the equation is a single number. Plugging into the spin-adapted formulas with the only allowed transition index $(i,a)=(g,u)$:

Two integrals appear and they mean different things physically:

• $(gu|gu)$ is the Coulomb self-repulsion of the transition density ${\psi }_g{\psi }_u$. This is the term that fires when the excitation creates a non-zero charge distribution that interacts with itself.
• $(gg|uu)$ is the Coulomb interaction between the occupied density ${\psi }_g^2$ and the virtual density ${\psi }_u^2$. Physically this is the electron-hole attraction — the exciton binding — that pulls the excitation energy below the bare orbital gap.

Tamm-Dancoff (TDA)

Drop $B$ and the singlet/triplet excitation energies are just the $A$ values:

The singlet-triplet splitting at the TDA level is therefore:

Twice the Coulomb self-repulsion of the transition density. That's Hund's rule made explicit in the algebra: the triplet pays no Coulomb cost for the excitation because its two electrons have parallel spins and the Coulomb interaction between same-spin transition density components cancels in the spin-adapted combination. The singlet doesn't get that cancellation and pays the full $2(gu|gu)$ penalty. Triplet wins, singlet loses.

Full Casida (TDHF)

Including $B$ and using the Hermitian rewrite ${\omega }^2=(A-B)(A+B)$ from the section above:

Two observations. First, both expressions sit at or below their TDA counterparts — the de-excitation coupling $B$ always lowers positive-frequency roots (this is the perturbative argument from the Try First probe below, now visible directly in the algebra). Second, the triplet expression's first factor goes negative when $(gu|gu)+(gg|uu)>({\epsilon }_u-{\epsilon }_g)$ — at that point ${\omega }_T^2$ goes negative, ${\omega }_T$ becomes imaginary, and you have a triplet instability. In the minimal basis you can see it happen as an inequality between three numbers.

From here to 6-31G

Everything below — the 6-31G PySCF demonstration, the singlet-triplet splitting, the TDA overshoot — is this same calculation done with three virtual orbitals instead of one. The Casida block becomes 3 × 3 instead of 1 × 1, and there's no closed form. But the structure is identical: diagonal of $A$ = orbital gap, plus Coulomb of the transition density, minus electron-hole attraction. The minimal-basis case is just every term reduced to a single integral.

H₂ in 6-31G: working through it by hand

H2 in a minimal-ish basis is the perfect testbed: one occupied orbital (${\sigma }_g$), a handful of virtuals, the whole Casida problem is a 3×3 in each spin block. Below: do an RHF SCF, build the four Casida tensors directly from the molecular ERIs, diagonalize, and check against PySCF's library TDHF. The two implementations should agree to machine precision because they are literally the same equation; the point of writing it out is that the equation is no longer a black box.

"""
Casida's equation for H2 (TDHF flavor).

We build the linear-response problem from scratch from a restricted
Hartree-Fock reference: form the Casida A and B matrices in the
occupied-virtual MO basis, project into the singlet and triplet
spin-adapted blocks, solve the pseudo-eigenvalue problem, and verify
against PySCF's built-in TDHF.

  [ A  B ][X]       [ 1   0 ][X]
  [ B  A ][Y] = w   [ 0  -1 ][Y]

When A - B is positive-definite (true at the HF minimum for stable
closed-shell molecules), this is equivalent to the Hermitian problem

  (A - B)^{1/2} (A + B) (A - B)^{1/2} Z = w^2 Z.

The Tamm-Dancoff approximation drops B and diagonalizes A directly.
TDDFT differs only in the two-electron kernel: replace -K_exchange by
the XC kernel f_xc. The matrix structure is identical, so the code is
identical down to the slot where you fill in the kernel.
"""

import numpy as np
from pyscf import gto, scf, tdscf, ao2mo

# ---- 1. Reference: restricted HF on H2 at R = 0.74 A --------------------
mol = gto.M(atom='H 0 0 0; H 0 0 0.74', basis='6-31g', verbose=0)
mf  = scf.RHF(mol).run()

nocc = mol.nelectron // 2
nmo  = mf.mo_coeff.shape[1]
nvir = nmo - nocc
eps  = mf.mo_energy
print(f"H2 / 6-31G  nocc = {nocc}, nvir = {nvir},  E_HF = {mf.e_tot:.6f} Ha")
print(f"HOMO-LUMO gap = {eps[nocc] - eps[nocc-1]:.6f} Ha "
      f"({(eps[nocc] - eps[nocc-1]) * 27.211386:.3f} eV)")

# ---- 2. Two-electron integrals (chemist notation) in the MO basis --------
eri_ao = mol.intor('int2e')
eri_mo = ao2mo.incore.full(eri_ao, mf.mo_coeff)

o = slice(0, nocc)
v = slice(nocc, nmo)

# All four-index tensors are stored with index order (i, a, j, b).
iajb = eri_mo[o, v, o, v]
ijab = eri_mo[o, o, v, v].transpose(0, 2, 1, 3)
ibja = eri_mo[o, v, o, v].transpose(0, 3, 2, 1)

# Orbital-energy gap on the diagonal of A: (eps_a - eps_i) delta_ij delta_ab
delta = np.zeros((nocc, nvir, nocc, nvir))
for i in range(nocc):
    for a in range(nvir):
        delta[i, a, i, a] = eps[nocc + a] - eps[i]

# ---- 3. Spin-adapted Casida A and B --------------------------------------
A_singlet = delta + 2.0 * iajb - ijab
B_singlet =         2.0 * iajb - ibja
A_triplet = delta              - ijab
B_triplet =                    - ibja

def casida_omega(A, B):
    """Solve [[A,B],[B,A]] [X;Y] = w [[1,0],[0,-1]] [X;Y].
    Returns the sorted positive excitation energies."""
    n = A.shape[0] * A.shape[1]
    A2 = 0.5 * (A.reshape(n, n) + A.reshape(n, n).T)
    B2 = 0.5 * (B.reshape(n, n) + B.reshape(n, n).T)
    ApB = A2 + B2
    AmB = A2 - B2
    w, U = np.linalg.eigh(AmB)
    sqrt_AmB = (U * np.sqrt(np.maximum(w, 0))) @ U.T
    M = sqrt_AmB @ ApB @ sqrt_AmB
    M = 0.5 * (M + M.T)
    omega2 = np.linalg.eigvalsh(M)
    return np.sort(np.sqrt(np.maximum(omega2, 0)))

def tda_omega(A):
    """Tamm-Dancoff: just diagonalize A."""
    n = A.shape[0] * A.shape[1]
    A2 = 0.5 * (A.reshape(n, n) + A.reshape(n, n).T)
    return np.sort(np.linalg.eigvalsh(A2))

om_S = casida_omega(A_singlet, B_singlet)
om_T = casida_omega(A_triplet, B_triplet)
om_S_tda = tda_omega(A_singlet)
om_T_tda = tda_omega(A_triplet)

# ---- 4. PySCF reference (full TDHF) --------------------------------------
nstates = nocc * nvir
tdhf_S = tdscf.TDHF(mf); tdhf_S.singlet = True;  tdhf_S.nstates = nstates; tdhf_S.kernel()
tdhf_T = tdscf.TDHF(mf); tdhf_T.singlet = False; tdhf_T.nstates = nstates; tdhf_T.kernel()
ref_S = np.sort(tdhf_S.e)
ref_T = np.sort(tdhf_T.e)

print(f"\nMax abs error vs PySCF TDHF: singlet {np.max(np.abs(om_S - ref_S)):.2e}, "
      f"triplet {np.max(np.abs(om_T - ref_T)):.2e}")

# ---- 5. Singlet-triplet splitting and TDA error --------------------------
Ha_eV = 27.211386
print("\nLowest excitations of H2 / 6-31G:")
print(f"  S0 -> T1   {om_T[0]:.4f} Ha = {om_T[0]*Ha_eV:6.3f} eV   (TDA: {om_T_tda[0]*Ha_eV:6.3f} eV)")
print(f"  S0 -> S1   {om_S[0]:.4f} Ha = {om_S[0]*Ha_eV:6.3f} eV   (TDA: {om_S_tda[0]*Ha_eV:6.3f} eV)")
print(f"  singlet - triplet splitting: {(om_S[0]-om_T[0])*Ha_eV:.3f} eV")
print(f"  TDA overshoots S1 by {(om_S_tda[0]-om_S[0])*Ha_eV:.3f} eV "
      f"(de-excitation coupling that TDA discards)")

Output:

H2 / 6-31G  nocc = 1, nvir = 3,  E_HF = -1.126755 Ha
HOMO-LUMO gap = 0.834290 Ha (22.702 eV)

Max abs error vs PySCF TDHF: singlet 2.89e-15, triplet 2.22e-15

Lowest excitations of H2 / 6-31G:
  S0 -> T1   0.3599 Ha =  9.793 eV   (TDA: 10.316 eV)
  S0 -> S1   0.5520 Ha = 15.020 eV   (TDA: 15.248 eV)
  singlet - triplet splitting: 5.226 eV
  TDA overshoots S1 by 0.228 eV (de-excitation coupling that TDA discards)

Agreement with PySCF is at the level of $10^{-15}$ Hartree across all three roots in each spin block — floating-point noise, as it should be. The triplet sits 5.2 eV below the singlet (Hund's rule in action), and the TDA bumps the singlet up by 0.23 eV.

Two things to flag about these numbers. First, they are HF-quality — TDHF, not TDDFT. The experimental H2 S0→T1 excitation is about 10.6 eV and S0→S1 (B¹Σu+) is about 12.6 eV; we're overshooting badly because HF ignores correlation and 6-31G has no diffuse functions for the diffuse excited states. Switching to TDDFT with a hybrid functional and a basis with diffuse augmentation (aug-cc-pVDZ or larger) brings it into experimental range. Second, this is the entire spectrum at this level — three singlet roots and three triplet roots, exhausting the $N_{\text{occ}}{N}_{\text{vir}}=3$ dimension of the problem. There is no information past that without enlarging the basis.

Why it doesn't scale to solids

For a molecule with a few hundred basis functions you get an $N_{\text{occ}}{N}_{\text{vir}}$ Casida matrix of size $\sim 10^4$, which dense LAPACK handles in seconds. For a solid in a plane-wave basis with a dense k-mesh the dimension blows up to $10^6$ and beyond, and worse, the spectrum becomes a continuum — there is no point resolving a million individual roots when what you want is the smooth function $\text{Im}\chi (\omega )$. That is the regime where Liouville-Lanczos takes over: never form the Liouvillian, never diagonalize it, get the whole spectrum from one matrix-free Krylov chain plus a continued-fraction terminator.

Deep dive: where A and B come from

Before diving into the algebra, an aside on the single-particle business that grinds my gears. Every single-particle model — Hartree-Fock, Kohn-Sham, the nuclear shell model — works by calculating a stack of single-particle levels and then filling them using some occupation rule (RHF, UHF, even shell, odd shell). You are approximating a many-body Hamiltonian by populating a single-particle level diagram. Single particle + spin rules is the formula. The fascinating part is that as soon as you get into "many-body" quantum mechanics you start learning how to make single-particle models in mean fields with spin rules.

Step 0: Strategy — linear response of the ground state

First, like most DFT calculations, you need the ground state. The DFT calculation returns with a set of orbitals which you fill from lowest to highest. In practice that means that the first $n/2$ orbitals are filled with two electrons each. Now we need to annoy, or as they say in the physics community perturb, the ground state. You know how when you annoy someone with just the right frequency they react? The same principle applies here. We use an external field and from there we test which frequencies annoy the ground state the most. Each frequency that gets a reaction is an excitation energy ${\omega }_n$, and the size of the reaction is the oscillator strength $f_n$.

Time to put that in math. The Greek letter $\delta$ just means "a small change in," so $\delta v_{\mathrm{ext}}$ is a small change in the external potential — that's our laser. We make it oscillate at frequency $\omega$:

The "c.c." at the end is shorthand for "complex conjugate of the previous term." We need it because the field has to be real, but $e^{-i\omega t}$ on its own is complex; adding its conjugate cleans that up. The electron density reacts to the field by wobbling — call that wobble $\delta \rho$, oscillating at the same frequency $\omega$.

How big a wobble depends on the frequency. To turn "how big" into a number, we use the dynamic polarizability $\alpha (\omega )$ — basically, how much the molecule's dipole moment shifts when you shake it with a field at frequency $\omega$:

Why bother writing it this way? Because $\alpha (\omega )$ blows up — has poles — at exactly the excitation energies. So if we can write down the equation that governs $\delta \rho (\omega )$, the eigenvalues of that equation are the ${\omega }_n$. That is the entire strategy.

Step 1: The TDKS equation

The starting point is the time-dependent Kohn-Sham equation for the one-particle density matrix:

Step 2: H_KS under a small perturbation

The Hamiltonian on the right is the instantaneous one, and it itself depends on the density. To first order in the perturbation, it splits into three pieces:

The first piece is the static KS Hamiltonian — its eigenstates are the ${\psi }_i$ and ${\psi }_a$, eigenvalues ${\epsilon }_i$ and ${\epsilon }_a$. The second is the external perturbation (the laser). The third is the self-consistent response: when the density wobbles by $\delta \rho$, the Hartree-XC potential wobbles with it, and the wobble acts back on the orbitals. The kernel $K_{\mathrm{HXC}}$ is the static response function — Coulomb minus exchange in TDHF, Coulomb plus the XC kernel $f_{xc}$ in TDDFT.

Step 3: Linear response of the density matrix → X and Y

Now expand the density matrix to first order: $\hat{\gamma }(t)={\hat{\gamma }}_0+{\hat{\gamma }}_1(t)$, where ${\hat{\gamma }}_0$ is the ground-state projector and ${\hat{\gamma }}_1$ is the wobble. The wobble has to live in the off-diagonal occupied-virtual block of the density matrix (the diagonal blocks can't change at first order — particle number is fixed). There are two independent corners: virtual-occupied (call its amplitude $X$) and occupied-virtual (amplitude $Y$):

Step 4: Insert into TDKS, go to frequency domain, project onto orbital pairs

Plug ${\hat{\gamma }}_1$ into the TDKS commutator and keep only first-order terms. The static-Hamiltonian piece gives the orbital-energy difference:

The self-consistent piece $[{\hat{V}}_{\mathrm{HXC}}^{\prime }[{\hat{\gamma }}_1],{\hat{\gamma }}_0]$ contributes the two-electron matrix elements — the kernel $K_{\mathrm{HXC}}$ sandwiched between orbital products. For a monochromatic field at frequency $\omega$, equate coefficients of $e^{-i\omega t}$ on both sides of $i{\partial }_t{\hat{\gamma }}_1=[\cdot ,\cdot ]$ and project onto orbital pairs $({\psi }_a,{\psi }_j)$. What drops out is exactly the Casida block equation, with matrix elements:

The bracket notation $⟨{\psi }_a{\psi }_j|K_{\mathrm{HXC}}|{\psi }_i{\psi }_b⟩$ means the spatial integral $\int {\psi }_a^{\ast }(\mathbf{r}){\psi }_j^{\ast }({\mathbf{r}}^{\prime })K_{\mathrm{HXC}}(\mathbf{r},{\mathbf{r}}^{\prime }){\psi }_i(\mathbf{r}){\psi }_b({\mathbf{r}}^{\prime })d\mathbf{r}d{\mathbf{r}}^{\prime }$ — the kernel acting on a product of two transition densities. That's the $(ia\Vert jb)$ notation used above. The difference between $A$ and $B$ is which orbital pair you put on each side of the kernel: $A$ pairs $(a,j)$ on the left with $(i,b)$ on the right; $B$ pairs $(a,b)$ with $(i,j)$. That's the only difference between excitation-excitation coupling and excitation-de-excitation coupling.

Step 5: TDDFT and TDHF — a structural difference in the coupling matrix

The derivation above applies to both TDDFT and TDHF — the only difference is the kernel $K_{\mathrm{HXC}}$. But the two cases have a hidden structural asymmetry worth surfacing. Because DFT's self-consistent potential is local (acts on the density at a single point), the TDDFT coupling matrix has the full set of orbital-index symmetries:

TDHF, whose exchange operator is nonlocal (it mixes the orbital at different spatial points), only enjoys the weaker hermiticity:

The practical consequence: adiabatic TDDFT matrices have fewer independent entries than TDHF ones, so they're cheaper to assemble and faster to diagonalize. Counterintuitive — you might expect TDHF to be simpler because it has no XC kernel to worry about — but locality buys symmetry, and symmetry buys speed. The same locality property is what makes the auxiliary-function (resolution-of-identity / density-fitting) factorization Casida introduced in his 1995 paper possible for TDDFT but harder for TDHF.

Step 6: Hermitian rewrite (covered in the body)

Having the block equation in hand, the next move is to fold its indefinite metric into the operator and end up with a standard Hermitian eigenvalue problem half the size: $(A-B{)}^{1/2}(A+B)(A-B{)}^{1/2}Z={\omega }^{2}Z$. The derivation — adding and subtracting the two block rows, using positive-definiteness of $A-B$ at a stable minimum to define the symmetric square root, then conjugating — is written out in the Why it's a pseudo-eigenvalue problem section above. No need to repeat it here.

Step 7: Spin adaptation (covered in the body)

For a closed-shell reference the full Casida problem factors into independent singlet and triplet blocks of half the size. The block-by-block A and B formulas — singlet picks up the $2(ia|jb)$ Coulomb piece that the triplet drops — appear in the Spin adaptation section above. The minimal-basis H₂ walk-through then shows the singlet-triplet splitting collapsing to a single integral.

Step 8: From derivation to code

Every step above has a one-to-one image in the H₂ Python program. Hovering the annotations on the code block lays this out line-by-line; the summary view:

• Steps 0–1 (reference state): scf.RHF(mol).run() produces the converged ${\psi }_i$, ${\epsilon }_i$, and density matrix. Everything downstream feeds off these.
• Step 2 (the kernel $K_{\mathrm{HXC}}$): In TDHF the kernel is the Coulomb-minus-exchange piece — computed from the two-electron integrals. mol.intor('int2e') + ao2mo.incore.full(eri_ao, mf.mo_coeff) assemble those integrals in the MO basis.
• Step 3 (X, Y amplitudes): Don't appear explicitly as variables in the program — they're implicit in the structure of the eigenvalue problem. The Casida matrix acts on a vector of length $N_{\mathrm{occ}}{N}_{\mathrm{vir}}$ which represents the combined X-Y degree of freedom in the spin-adapted block.
• Step 4 (assembling A and B): The four lines A_singlet = delta + 2.0 * iajb - ijab, B_singlet = 2.0 * iajb - ibja, A_triplet = delta - ijab, B_triplet = -ibja are this step, written out for both spin sectors.
• Step 5 (coupling-matrix symmetry): Doesn't appear directly in the code — it's why we get away with storing iajb, ijab, ibja as three small tensors rather than enumerating all asymmetric variants.
• Step 6 (Hermitian rewrite): The body of casida_omega. np.linalg.eigh(AmB) diagonalizes $A-B$; (U * np.sqrt(np.maximum(w, 0))) @ U.T builds the symmetric square root; sqrt_AmB @ ApB @ sqrt_AmB conjugates; np.linalg.eigvalsh(M) gives ${\omega }^2$; the final np.sqrt recovers $\omega$.
• Step 7 (spin adaptation): The reason there are two calls to casida_omega, one for each of (A_singlet, B_singlet) and (A_triplet, B_triplet). The blocks don't talk to each other; we just solve each independently.

Read the code once with this map in mind and the whole program is just the derivation, executed.

Q&A

Study notes in question-and-answer form. The questions are the ones a reader would actually ask at each point — framing first, then drill-down, then a final compression. Pair them with the deep dive above when you want to test recall.

What problem does Casida solve, and why was it hard before?

Computing electronic excitation energies and oscillator strengths. The hard way is to build excited-state many-body wavefunctions directly (CIS, CISD, CASSCF, coupled cluster) and compute transitions between them. Those methods scale steeply — $O(N^4)$ at the cheapest, $O(N^7)$ for chemically accurate ones. Casida sidesteps it: instead of building excited states, perturb the ground state with a small external field and look at the response. The poles of the response function are the excitation energies, and they come out of a single eigenvalue problem of dimension $N_{\mathrm{occ}}\cdot N_{\mathrm{virt}}$ — much smaller than the full many-body Hilbert space.

What's the core idea?

Linearize the time-dependent Kohn-Sham (or Hartree-Fock) equation around the converged ground state. The first-order density-matrix perturbation ${\hat{\gamma }}_1$ has only off-diagonal entries in the occupied-virtual basis — two amplitudes per particle-hole pair, called $X$ and $Y$. In frequency domain they obey a coupled block equation:

Find the eigenvalues $\omega$; those are the excitation energies.

Why does the dimension come out to N_occ × N_virt?

$X$ and $Y$ are indexed by particle-hole pairs $(i,a)$ — $i$ over occupied orbitals, $a$ over virtual orbitals. There are $N_{\mathrm{occ}}\cdot N_{\mathrm{virt}}$ such pairs. For H₂ in 6-31G with one occupied and three virtuals, that's 3. For a medium molecule in cc-pVDZ with 100 occupied and 200 virtual, that's 20,000 — solvable by dense LAPACK in seconds. For a solid with thousands of bands and a dense k-mesh the dimension explodes past $10^6$. That's where Casida breaks down and Liouville-Lanczos takes over.

What are X and Y physically?

$X_{ia}$ is the amplitude for an electron to leave occupied orbital $i$ and end up in virtual orbital $a$ at frequency $\omega$. $Y_{ia}$ is the amplitude for the time-reversed process — an electron returning from $a$ to $i$. Linear response includes both because the density's wobble at $\omega$ has a forward and a backward component. They live in the two off-diagonal corners of ${\hat{\gamma }}_1$; the diagonal blocks are forbidden by particle-number conservation and idempotency.

Why does the block equation have a minus sign on the right?

$X$ and $Y$ are time-reversed partners. An excitation at $+\omega$ is the same physical state as a de-excitation at $-\omega$. The minus sign in the right-hand metric tracks this. Positive eigenvalues are physical excitations; negative eigenvalues are their de-excitation partners (same states, opposite sign of $\omega$). The indefinite metric is what keeps Casida from being a standard eigenvalue problem until you do the Hermitian rewrite.

How does the Hermitian rewrite work?

Add the two block rows: $(A+B)(X+Y)=\omega (X-Y)$. Subtract: $(A-B)(X-Y)=\omega (X+Y)$. Combine: $(A+B)(A-B)(X-Y)={\omega }^2(X-Y)$. The operator $(A+B)(A-B)$ isn't Hermitian, but if $A-B$ is positive-definite (which it is at a stable minimum), take its symmetric square root and conjugate:

Hermitian, positive-definite, half the size of the original. Eigenvalues are ${\omega }^2$; physical excitation energies are the square roots.

What's the difference between TDA and full Casida?

TDA drops $B$ entirely. The Casida equation collapses to $AX=\omega X$ — a standard Hermitian eigenvalue problem of half the size, no rewrite needed. The cost is a systematic upward bias on positive-frequency roots, because adding $B$ always lowers them. The benefit: TDA is immune to triplet instabilities (no $A-B$ to go non-positive-definite), and the eigenvectors are orthonormal under the standard inner product. Most molecular excitation-energy work uses TDA by default; full Casida is what you reach for when oscillator strengths or specific frequency regions matter.

What does spin adaptation give you for closed-shell?

The full spin-orbital problem factors into independent singlet and triplet blocks, each half the size. Matrix elements:

The singlet picks up the $2(ia|jb)$ Coulomb-of-transition-density piece that the triplet drops. Both keep the orbital gap and the electron-hole attraction $-(ij|ab)$. Solve each block independently; they don't mix.

Why is the singlet-triplet splitting exactly 2(gu|gu) in minimal-basis H₂?

Minimal basis gives one occupied (${\psi }_g$) and one virtual (${\psi }_u$). The only allowed index is $(g,u)$, so the spin-adapted $A$ matrices are single numbers: $A^S=({\epsilon }_u-{\epsilon }_g)+2(gu|gu)-(gg|uu)$, $A^T=({\epsilon }_u-{\epsilon }_g)-(gg|uu)$. Difference: $2(gu|gu)$. Physically: the singlet pays a Coulomb cost the triplet avoids by aligning spins (exchange hole). Hund's rule, written out in algebra.

What's a triplet instability?

A point where the lowest eigenvalue of $A-B$ in the triplet block crosses zero. Past it, ${\omega }^2<0$ and the lowest triplet excitation energy goes imaginary. Physically: the closed-shell RHF reference becomes unstable to spin-symmetry breaking — a small density wobble in that direction lowers the total energy. Classic example: H₂ stretched past the Coulson-Fischer point, where the RHF reference no longer correctly describes the dissociating bond. The fix is to relax the closed-shell constraint and use a UHF (or broken-symmetry) reference.

How do oscillator strengths come out of the eigenvectors?

After the Hermitian rewrite gives eigenvalues ${\omega }_n^2$ and eigenvectors $Z_n$, undo the rewrite to recover $(X+Y{)}_n\propto (A-B{)}^{-1/2}{Z}_n$. The transition density is ${\rho }_n(\mathbf{r})={\sum }_{ia}(X+Y{)}_{n,ia}{\psi }_i^{\ast }(\mathbf{r}){\psi }_a(\mathbf{r})$. The oscillator strength is

Dipole-integrate the transition density, square, multiply by $2{\omega }_n/3$. Dimensionless number that says how brightly state $n$ absorbs.

Why does pure DFT-TDDFT badly underestimate charge-transfer energies?

At long range, the xc potential $v_{xc}$ of LDA and GGA decays exponentially instead of the correct $-1/r$. For a charge-transfer excitation between well-separated donor and acceptor, the right asymptote is $I_D-A_A-1/R$: donor ionization potential, minus acceptor electron affinity, minus the Coulomb attraction of the resulting ion pair. Adiabatic TDDFT collapses the CT energy to roughly ${\epsilon }_a-{\epsilon }_i$, which has no $1/R$ piece — the missing long-range exchange was supposed to provide it. HF exchange is nonlocal and has the correct $-1/r$ tail by construction, which is why TDHF gets CT right and why range-separated hybrids (CAM-B3LYP, $\omega$B97X, LC-$\omega$PBE) restore long-range exchange to fix this.

Why does Casida fail for solids?

Two reasons. First, dimension: in a plane-wave calculation with thousands of bands and a dense k-mesh, the Casida matrix dimension easily reaches $10^6$ or more. Dense diagonalization is impossible at that scale. Second, the spectrum: a solid has a continuum, not a discrete set of poles. Resolving each eigenvalue individually is meaningless — what you want is the smooth function $\mathrm{Im}\chi (\omega )$. Liouville-Lanczos handles both: it computes the response function directly via a matrix-free Krylov chain, without ever forming or diagonalizing the Liouvillian.

What's the cleanest mental compression?

Two steps. (1) Linear response: perturb the converged ground state with a small periodic field; the density wobbles; the poles of the response function are the excitation energies. (2) Eigenvalue problem: project the linearized TDKS equation onto particle-hole orbital pairs; get a block matrix equation with an indefinite metric; Hermitian-rewrite it to size $N_{\mathrm{occ}}{N}_{\mathrm{virt}}$; eigenvalues are ${\omega }^2$. Everything else — spin adaptation, TDA, oscillator strengths, triplet instabilities, charge-transfer failure, the auxiliary-function trick — is a consequence of those two steps.

Related on this site

Liouville-Lanczos is the matrix-free, frequency-sweep cousin of Casida for solids and large systems. Lanczos iteration is the underlying Krylov method; bi-orthogonal Lanczos is the non-Hermitian variant needed when the operator (the Liouvillian) is non-symmetric.

Unit conversion used throughout: $1\mathrm{Ha}=27.211386\mathrm{eV}$ (CODATA 2018 recommended value).