Transforming AO basis electron repulsion integrals into molecular spin orbital basis?

Question

I am attempting to do project 5 on Daniel Crawford's website. I'm having some trouble reproducing the MP2 energy, which is given in project 4.

In project 4, the MP2 energy is given as:

$$ E_{MP2} = \sum_{ij}\sum_{ab} \frac{(ia|jb)[2(ia|jb)-(ib|ja)]}{e_i+e_j - e_a - e_b}$$

Where i, j are occupied orbitals and a, b are unoccupied orbitals. The ERIs given above are in the MO basis. Crawford's website provides the ERIs in the AO basis. It describes that you can convert the AO basis ERIs into the MO basis ERIs by summing over the MO C coefficients, which I have prepared with SCF, in this fashion:

$$ (pq|rs) = \sum_{\mu\nu\lambda\sigma} C^p_\mu C^q_\nu C^r_\lambda C^s_\sigma (\mu\nu|\lambda\sigma) $$

Where $p, q, r, s$ are MO basis ERIs and $\mu \nu \lambda \sigma $ are AO basis ERIs. (This is the inefficient formulation of this sum but it looked easier to type in LaTeX; the better version is on the project 4 page which is what I actually coded). I was able to prepare the MO basis ERIs and the MP2 energy correctly (the website has the right answer) in this method.

On to project 5, where my problem is. Project 5 is to make a CCSD code which uses the spin MO basis as opposed to the non-spin MO basis that MP2 uses, but an intermediate part of that is to show that you can calculate the MP2 energy using these spin MO orbitals. However, I am having trouble doing so. I suspect my problem is that I am not converting my ERIs properly, but I am not 100% sure.

I am trying to use this formula for the MP2 energy (from project 5): $$E_{MP2} = \frac{1}{4}\sum_{ijab}\left<ij||ab\right> t^{ab}_{ij}$$

I believe it is the same as project 4, where $i,j$ are occupied orbitals, $a,b$ are unoccupied. The formula for $t^{ab}_{ij} $ is given on the project 5 page. When I calculate the MP2 energy with the above formula I get something quite different than what I get from the first MP2 energy formula (from project 4)

I think my problem is that I am misinterpreting what $\left<ij||ab\right>$ means, but I am not sure in what way. The site provides an example code to prepare these (given here), but it's in C++ and I am just doin this in python to practice, so I am making my own. Also, following the advice on the page, I am storing all the ERIs (AO basis, MO basis, spin MO basis) as large 4 index tensors and not using any fancy indexing, just to make it easier to debug. The example code seems simple enough that I think I have reproduced it correctly (except for maybe storing the results), so I think my problem is not a coding one, rather that I am misunderstanding some idea or notation. I am doing three things to obtain it from the given AO ERIs:

1. Construct the MO ERIs from the AO ERIs, identically to in project 4 (see the four index sum over greek letters I put above). I know I am doing this right as I used these to get the MP2 energy in project 4 correctly.

2. Construct the spin MO basis ERIs. I use project 5's advice of treating spin as the number of the function (ordered by energy). I have: $$ \left<pq|rs\right> = (pr|qs) \delta_{r\%2, p\%2} \delta_{s\%2, q\%2} $$ Where $\left<pq|rs\right>$ is the spin MO integral and $(pr|qs)$ is the spatial MO integral. I have the orbitals sorted by energy (the energy given by SCF) and numbered. Say p is #2 and r is #3; then p mod 2 = 0 and r mod 2 = 1, so p would be spin up and r would be spin down, so the spin orbital is killed. This is the method that project 5 is using to simulate spin (at least, that's how I'm interpreting the example code. Could that be my error?).

3. I anti-symmetrize the spin MO integrals like this: $$ \left<pq||rs\right> = (pr|qs) - (ps|qr) $$ I didn't write them here, but the two integrals on the right side have the spin deltas built in to my code as well. Since all the ERIs in the CCSD equations are anti symmetrized ones, I store the 4 index ERI tensor in this form.

I think my problem may be that I am unconsciously mixing up physicists and chemists notation somewhere in my code when I do step 2 in my conversion? In project 5's example code, they use an indexing scheme to read in the spatial MO ERIs, which I have not implemented or really understood, but I think I am emulating correctly by just using a big 4 index tensor to store the ERIs. In my code, I am reading in $(pr|qs)$ and $(ps|qr)$ like this:

prqs = ERIs[p][r][q][s] * (p%2 == r%2) ...
psqr = ERIs[p][s][q][r] 

I didn't type all of them but the p%2 == r%2 thing is meant to kill it if the spins are not equal (aka, kill it if the number of the index of the orbitals are of different parity). Here, prqs is in chemists notation (which is what I assume they mean by using parenthesis) and ERIs is my four index tensor (using numpy) of the spatial MO ERIs. Then I antisymmerize them and store them like this:

spin_ERIs[p][q][r][s] = prqs - psqr

by spin_ERIs[p][q][r][s], I mean $\left<pq||rs\right>$. Could my problem be that since I am switching between chemists and physicists notation (again, I assume the <> means physicists notation and () chemists) I need to switch some indices around, or something? It seems like that is built in already, since we have $\left<pq||rs\right> = (pr|qs) - (ps|qr)$, not $\left<pq||rs\right> = \left<pq|rs\right> - \left<qp|rs\right> $or something, but I'm not sure if that's the problem, or something else entirely.

I guess this problem is pretty convoluted but basically, can somebody tell me that I am thinking the right way from my steps 1, 2, 3? Thanks.

Answer 1

I am not sure whether this is it, but I think your equation in step 2 should be:

\begin{equation} \langle pq|rs\rangle = \left(\frac{p}{2}\frac{r}{2}\middle|\frac{q}{2}\frac{s}{2}\right) \delta_{r\%2, p\%2} \delta_{s\%2, q\%2} \end{equation}

The idea for $(pr|qs)$ is to store $\frac{N}{2}$ spatial orbitals representing $N$ spin orbitals, by remembering that each orbital appears twice (one $\alpha$, one $\beta$). So the indices $p,q,r$ and $s$ go only up to $\frac{N}{2}$ for the spatial ERIs.

Note that we use an integer division here: $\frac{3}{2}=1$ and $\frac{2}{2}=1$. Thus there are always two spin indices mapping to the same spatial index. One might also write the spin indices as compound indices, consisting of one spatial index and the spin state: $0\alpha, 0\beta, 1\alpha, 1\beta, \dots$.

But for the spin ERIs $\langle pq|rs\rangle$ the indices run up to $N$. So with $(pr|qs)$ you would try to access elements outside of it (as soon as the indices become greater then $\frac{N}{2}$). You can also find this in the code of your linked hint.

Furthermore, I think you are mixing up the different notations. Hopefully this will clarify it: According to the book by Szabo and Ostlund the notation for these integrals should be:

• $\langle ij|kl\rangle$ is physicists' notation for spin orbitals ($i$ and $j$ are bra-vectors of different electronic coordinates)
• $[ij|kl]$ is chemists' notation for spin orbitals ($i$ and $j$ are bra- and ket vector of the same electronic coordinate)
• $(ij|kl)$ denotes spatial orbitals following the ordering of the chemists' notation

There is no "physicists'-like" notation for spatial orbitals.

We can relate chemists' and physicists' notation by

\begin{equation} [ij|kl] = \left<ik|jl\right> \end{equation}

Also it defines \begin{equation} \langle pq||rs\rangle = \langle pq|rs\rangle - \langle qp|rs\rangle \end{equation}

which means your

\begin{equation} \langle pq||rs\rangle = (pr|qs) - (ps|qr) \end{equation}

should be wrong. But maybe you are already using the correct thing in your code there, just with the wrong notation?

Answer 2

I think the answer given for Fedoran is what you are asking for, on a site note, for the implementation you are doing take a look at Coupled Cluster with Singles and Doubles (CCSD) in Python , I think I found a few mistakes in it, but it gives the overall most simple implementation. A tested to work implementation that is very equivalent to the previous link can be found here, python CCSD(T)