11
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Have you tried constraining yourself to closed shell RHF only? It would eliminate the complications due to spin and if you get right results there, you know your problem is in your handling of the spin orbitals. $\endgroup$ – Fl.pf. Jul 7 '17 at 7:52
  • $\begingroup$ I am doing so. The test case on the website is water, using STO-3G. So 7 basis functions, 10 electrons. 5 lowest orbitals doubly filled, top 2 are unoccupied. $\endgroup$ – iammax Jul 7 '17 at 13:41
7
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Well, I used: $$ <pq||rs> = (pr|qs) - (ps|qr) $$ Because that seems to be what the "hint script" is indicating; that that I mean, that is what it uses for "value1" and "value2". I suppose I am not 100% sure that the original integrals provided (the spatial AO basis integrals) are in chemist's notation or not; I assume they are because they use (), but maybe that's just the notation on the site? (reply 1/2) $\endgroup$ – iammax Jul 7 '17 at 18:57
  • 1
    $\begingroup$ Yes, the spin orbitals would be something like "0a, 0b, 1a, 1b, ..." but since we can only index with integers this has to be translated to "0, 1, 2, 3, ...". It is important here that the $\frac{p}{2}$ division will be an integer division, e.g. 3/2 = 1 and 2/2=1, so there are always two spin indices mapping to the same spatial index. $\endgroup$ – Feodoran Jul 7 '17 at 19:08
  • 1
    $\begingroup$ About your first comment: You're right, the code in the hint already subtracts the exchange integral. Which means it takes care of step 2 AND step 3. $\endgroup$ – Feodoran Jul 7 '17 at 19:11
  • 2
    $\begingroup$ Your energy assignments are correct, but keep in mind that spatial orbitals can be doubly occupied and spin orbitals only singly. And we do not change the number of electrons in the system. So 10 electrons can fill 5 spatial orbitals, but we need 10 spin orbitals (0 to 9) for them! $\endgroup$ – Feodoran Jul 7 '17 at 20:38
  • 1
    $\begingroup$ Ah, and with that (remembering that the spin orbitals are 1-slot, not 2-slot), I am getting the right MP2 energy... can't believe I overlooked such obvious things, but thank you. $\endgroup$ – iammax Jul 7 '17 at 20:55
5
$\begingroup$

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)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.