What are the other matrix elements in singlet symmetry-adapted CISD?

Question

Szabo/Ostlund list the CI matrix elements between singlet symmetry-adapted configurations (SAC) in Table 4.1 of their book:

$$ \langle ^1\Psi^r_ a \lvert \mathcal{H} - E_ 0 \rvert ^1\Psi^r_ a \rangle, \langle \Psi_ 0 \lvert \mathcal{H} \rvert ^1\Psi^{rr}_ {aa} \rangle, \langle \Psi_ 0 \lvert \mathcal{H} \rvert ^1\Psi^{rs}_ {aa} \rangle, \langle \Psi_ 0 \lvert \mathcal{H} \rvert ^1\Psi^{rr}_ {ab} \rangle, \langle \Psi_ 0 \lvert \mathcal{H} \rvert ^A\Psi^{rs}_ {ab} \rangle, \langle \Psi_ 0 \lvert \mathcal{H} \rvert ^B\Psi^{rs}_ {ab} \rangle, \langle ^1\Psi^{rr}_ {aa} \lvert \mathcal{H} - E_ 0 \rvert ^1\Psi^{rr}_ {aa} \rangle, \langle ^1\Psi^{rs}_ {aa} \lvert \mathcal{H} - E_ 0 \rvert ^1\Psi^{rs}_ {aa} \rangle, \langle ^1\Psi^{rr}_ {ab} \lvert \mathcal{H} - E_ 0 \rvert ^1\Psi^{rr}_ {ab} \rangle, \langle ^A\Psi^{rs}_ {ab} \lvert \mathcal{H} - E_ 0 \rvert ^A\Psi^{rs}_ {ab} \rangle, \langle ^B\Psi^{rs}_ {ab} \lvert \mathcal{H} - E_ 0 \rvert ^B\Psi^{rs}_ {ab} \rangle, \langle ^A\Psi^{rs}_ {ab} \lvert \mathcal{H} \rvert ^B\Psi^{rs}_ {ab} \rangle. $$

The expressions are given in terms of restricted, canonical MOs. While I have been able to arrive at the same results for all listed elements I tried so far, I am apparently unable to correctly derive elements that are not listed. I do not believe that they are all $0$.

I am testing my toy SCF/MP2/CIS/CID implementation against ORCA and am able to reproduce RHF, UHF, RMP2, and RCIS/TDHF results for different systems with good precision. However, CID in general eludes me thus far. $\ce{H_2}$ in a minimal basis set (single-$\zeta$) works correctly, as does $\ce{He}$ in double-$\zeta$. However, my results for $\ce{He}$ in triple-$\zeta$ are rather far off.

I am looking for the correctly derived off-diagonal elements of the CISD matrix. Lacking this, concrete pointers on which other, freely available QC suite will print the matrix are also welcome. Unfortunately, I am reduced to the status of a hobbyist without access to literature.

Edit: In the comments, it was asked how I was confident in the integrals. The AO integral code is ancient and has been verified for RHF and UHF against ORCA and Turbomole. While performing the CID calculation, the RHF and RMP2 energies are being calculated at the same time, and match the ORCA results. I just performed a calculation for $\ce{BeH2}$ with matching results - given the number of AO/MO involved, I feel confident in my AO-MO-transformation.

---

A Szabo, NS Ostlund Modern Quantum Chemistry, Dover Publications, first edition, 1996.

Answer 1

An easier route to using symmetry-adapted CI with hand-derived matrix elements is to implement CI with determinant strings. That is, you construct bitstrings of which orbitals are occupied in the determinant, and you don't care about adapting your basis for $\hat{S}^2$. This is the way most codes work, since the resulting algorithm is easy to make very fast, whereas spin-adaptation may require you to limit your expansion length.

The string based CI works in the space of spin-orbitals. For instance, the Hartree-Fock state would be $|{\rm HF}\rangle = (1, 1, \dots, 1, 1, 0, 0, \dots, 0, 0)^{\rm T}$, and the first excited determinant would be $(1, 1, \dots, 1, 0, 1, 0, \dots, 0, 0)^{\rm T}$. What you need to do is just build all the determinants that have the wanted $\hat{S}_z$, i.e. the correct number of $\alpha$ and $\beta$ electrons.

Now, building the CI Hamiltonian is very straightforward: once you have the bitstrings $|i\rangle$, you can obtain the matrix elements $\langle i | H | j \rangle$ with e.g. the Slater-Condon rules. The matrix element is zero unless $|i\rangle$ and $|j\rangle$ differ by at most a double excitation; you can find out the difference in the states' occupations with a bit-wise XOR operation. Counting the bits in the result of the XOR gives 0 if the determinants match, 2 if they differ by a single excitation, 4 in the case of a double excitation, and you don't care about the rest since the matrix element is zero.

You can build the sparse Hamiltonian in memory, if you don't go to huge active spaces just by looping over the pairs of bit strings, and use library implementations of sparse matrix diagonalizers.

The best thing about using this method, in addition to being relatively straightforward to implement, is that you can do arbitrary level CI with it: CID, CISD, CISDT, CISDTQ, ..., all the way to FCI and CASSCF.