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

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

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$.

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.

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, $AX=\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\times 10^{-9}$ Ha, about what is left of the self-consistent field's convergence; the tests allow $10^{-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 =\left(\begin{array}{cc}A& B\\ B& A\end{array}\right)$, the lowest TDHF $\omega$ is the minimum of $(X,Y{)}^T\Lambda (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}AX/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|jb)$, whose diagonal $2(ia|ia)$ is the Coulomb self-repulsion of the charge distribution ${\varphi }_i{\varphi }_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.