# On the spin-adaptation of 2nd order Møller-Plesset perturbation theory

In spin adaptation (common in electronic structure theory) we approximate that the spin-up and spin-down spin orbitals (herein assumed to be real and represented by $$\chi$$ with some placeholder spin coordinate $$\vec{s}$$ and spatial coordinates $$\vec{x}$$) have the same spatial orbitals (represented by $$\phi$$ with spatial coordinates $$\vec{x}$$). This allows us to eliminate the spin components of the molecular orbitals (MOs) in exchange for an 8-fold speed-up of two-electron integral calculation and storage:
$$\chi(\vec{x_1})=\phi(\vec{r_1})*\omega(\vec{s_1}) \rightarrow \phi(\vec{r_1})$$
I am trying to spin-adapt the spin-orbital, second-order MP2 energy correction (using canonical MOs):
$$E^{(2)}=\frac{1}{4}\sum_{ijab}^{SO} \frac{\Bra{ij}\Ket{ab}^{2}}{\Delta _{ab}^{ij}}$$
(where $$\Delta _{ab}^{ij}=\epsilon{_i}+\epsilon{_j}-\epsilon{_a}-\epsilon{_b}$$, where i and j represent occupied MOs in the reference solution, where a and b represent unoccupied MOs in the reference solution, and where $$\epsilon_m$$ denotes the orbital energy of the $$m^{th}$$ molecular orbital) to arrive at the closed-shell (spin-adapted) expression:
$$E^{(2)}=\sum_{ijab}^{SF} \frac{2\Bra{ij}{ab}\left.\right>^{2}-\Bra{ij}{ba}\left.\right>\Bra{ij}{ab}\left.\right>}{\Delta _{ab}^{ij}}$$

Below is my attempted derivation (I use SO and SF to denote whether the indices i,j,a,b denote spin-orbital or spin-free MOs, respectively):

$$E^{(2)}=\frac{1}{4}\sum_{ijab}^{SO} \frac{\Bra{ij}\Ket{ab}^{2}}{\Delta _{ab}^{ij}}$$
$$=\frac{1}{4}\sum_{ijab}^{SF} \frac {\Bra{i_\alpha j_\alpha}\Ket{a_\alpha b_\alpha}^{2} +\Bra{i_\alpha j_\beta}\Ket{a_\alpha b_\beta}^{2} +\Bra{i_\alpha j_\beta}\Ket{a_\beta b_\alpha}^{2} +\Bra{i_\beta j_\alpha}\Ket{a_\alpha b_\beta}^{2} +\Bra{i_\beta j_\alpha}\Ket{a_\beta b_\alpha}^{2} +\Bra{i_\beta j_\beta}\Ket{a_\beta b_\beta}^{2} } {\Delta _{ab}^{ij}}$$
(note: there are other possible terms, but they equate to zero due to orthonormality of spin component. Also note that I have assumed spin-up and spin-down MO energies are equivalent, allowing me to write only one denominator)
We now expand the anti-symmetrized two electron integrals using $$\Bra{ij}\Ket{ab}=\Bra{ij}\left.ab\right> - \Bra{ij}\left.ba\right>$$ and the orthonormality of the spin component, $$\Bra{i_{\omega_{1}}}\left.i_{\omega_{2}}\right>=\delta_{\omega_{1},\omega_{2}}$$, keeping only non-zero terms:

$$E^{(2)}=\frac{1}{4}\sum_{ijab}^{SF} \frac {{\Bra{ij}\Ket{ab}}^{2} +\Bra{ij}\left.ab\right>^{2} +\Bra{ij}\left.ba\right>^{2} +\Bra{ij}\left.ba\right>^{2} +\Bra{ij}\left.ab\right>^{2} +{\Bra{ij}\Ket{ab}}^{2} } {\Delta _{ab}^{ij}}$$

Where I have kept the ordering of terms to facilitate term-by-term comparison to the previous energy expression. We can, of course, combine these terms:

$$E^{(2)}=\frac{1}{4}\sum_{ijab}^{SF} \frac {2{\Bra{ij}\Ket{ab}}^{2} +2\Bra{ij}\left.ab\right>^{2} +2\Bra{ij}\left.ba\right>^{2}} {\Delta _{ab}^{ij}} = \frac{1}{2}\sum_{ijab}^{SF} \frac {{\Bra{ij}\Ket{ab}}^{2} +\Bra{ij}\left.ab\right>^{2} +\Bra{ij}\left.ba\right>^{2}} {\Delta _{ab}^{ij}}$$

Almost there! Now we expand the anti-symmetrized two electron integrals to get:

$$E^{(2)}= \frac{1}{2}\sum_{ijab}^{SF} \frac {\Bra{ij}\left.ab\right>^{2} -2\Bra{ij}\left.ab\right>\Bra{ij}\left.ba\right> +\Bra{ij}\left.ba\right>^{2} +\Bra{ij}\left.ab\right>^{2} +\Bra{ij}\left.ba\right>^{2}} {\Delta _{ab}^{ij}}$$

Which, upon combining terms and factoring out a 2, yields:

$$E^{(2)}= \sum_{ijab}^{SF} \frac {\Bra{ij}\left.ab\right>^{2} -\Bra{ij}\left.ab\right>\Bra{ij}\left.ba\right> +\Bra{ij}\left.ba\right>^{2} } {\Delta _{ab}^{ij}}$$

This is SO close, but not equal, to the spin-adapted energy expression for MP2 (reproduced below):

$$E^{(2)}=\sum_{ijab}^{SF} \frac{2\Bra{ij}{ab}\left.\right>^{2}-\Bra{ij}{ba}\left.\right>\Bra{ij}{ab}\left.\right>}{\Delta _{ab}^{ij}}$$

But we're in trouble, as $$\Bra{ij}{ab}\left.\right>^{2} \neq \Bra{ij}{ba}\left.\right>^{2}$$, so the above expressions apparently do not equate. Additionally, two-electron integral symmetries do not help us out. At the time of writing, I had not found the solution. But in the nick of time, I found the workaround, and have posted it below.

It turns out the trick is not a matter of exploiting symmetry or alternate derivations, but rather a simple bookkeeping trick!

$$E^{(2)}= \sum_{ijab}^{SF} \frac {\Bra{ij}\left.ab\right>^{2} -\Bra{ij}\left.ab\right>\Bra{ij}\left.ba\right> +\Bra{ij}\left.ba\right>^{2} } {\Delta _{ab}^{ij}}$$ can be broken into three sums:

$$E^{(2)}= \sum_{ijab}^{SF} \frac {\Bra{ij}\left.ab\right>^{2}} {\Delta _{ab}^{ij}} -\sum_{ijab}^{SF} \frac{\Bra{ij}\left.ab\right>\Bra{ij}\left.ba\right> } {\Delta _{ab}^{ij}} +\sum_{ijab}^{SF} \frac{\Bra{ij}\left.ba\right>^{2} } {\Delta _{ab}^{ij}}$$

The beauty is that the indices of each summation are independent of the indices of the other summations. As such, I am free to rename a as b and b as a in the third summation:

$$\sum_{ijab}^{SF} \frac{\Bra{ij}\left.ba\right>^{2} } {\Delta _{ab}^{ij}}= \sum_{ijba}^{SF} \frac{\Bra{ij}\left.ab\right>^{2} } {\Delta _{ba}^{ij}}$$

Which, because a and b behave the same in the summation (both over unoccupied MOs) and in the denominator (both represent the energy of an unoccupied MO), can be rearranged as:

$$\sum_{ijab}^{SF} \frac{\Bra{ij}\left.ab\right>^{2} } {\Delta _{ab}^{ij}}$$

(In pseudo math, $$\sum_{ijba}^{SF} = \sum_{ijab}^{SF}$$ and $${\Delta _{ab}^{ij}}= {\Delta _{ba}^{ij}}$$). Putting it all together, we get

$$E^{(2)}= \sum_{ijab}^{SF} \frac {\Bra{ij}\left.ab\right>^{2}} {\Delta _{ab}^{ij}} -\sum_{ijab}^{SF} \frac{\Bra{ij}\left.ab\right>\Bra{ij}\left.ba\right> } {\Delta _{ab}^{ij}} +\sum_{ijab}^{SF} \frac{\Bra{ij}\left.ab\right>^{2} } {\Delta _{ab}^{ij}} = 2\sum_{ijab}^{SF} \frac {\Bra{ij}\left.ab\right>^{2}} {\Delta _{ab}^{ij}} -\sum_{ijab}^{SF} \frac{\Bra{ij}\left.ab\right>\Bra{ij}\left.ba\right> } {\Delta _{ab}^{ij}} = \sum_{ijab}^{SF} \frac{2\Bra{ij}{ab}\left.\right>^{2}-\Bra{ij}{ba}\left.\right>\Bra{ij}{ab}\left.\right>}{\Delta _{ab}^{ij}}$$

Which proves the equality.

Spin adaptation is certainly a valuable thing in correlation calculations. Unfortunately, I only know of simple schemes for 2 different cases: (1) closed-shell singlets ($$S = 0$$, $$m_S = 0$$) and (2) closed shell triplets ($$S = 1$$, $$m_S = 0$$). The formula you have cited is closed-shell singlet adaptation. Since the above answers prove the result algebraically, I thought I could show the same result using diagrams which are a very useful tool for, among other things, deriving MPn energy and wavefunction corrections. The Hamiltonian used in post-Hartree Fock correlation calculations is the normally-ordered form

$$H_N = F_N + V_N = \sum_{p,q = 1}^{\#\,\mathrm{spinorbitals}} \langle p | f | q \rangle N[a^p a_q] + \dfrac{1}{2} \sum_{p,q,r,s = 1}^{\#\,\mathrm{spinorbitals}} \langle pq | rs \rangle N[a^p a^q a_s a_r]$$

where $$a^p_{(p)}$$ are the fermionic creation (annhiliation) operators, $$H_N = H - \langle \Phi_0 | H | \Phi_0 \rangle$$, $$|\Phi_0\rangle$$ is the Hartree-Fock reference Slater determinant, $$N[\ldots]$$ denotes normal ordering with respect to the Fermi vacuum $$|\Phi_0\rangle$$, and $$\langle p | f | q \rangle$$ and $$\langle pq | rs \rangle$$ denote the usual Fock and 2-electron integral matrix elements in the MO basis. The expectation value $$\langle \Phi_0 | H | \Phi_0 \rangle$$ is simply the Hartree-Fock electronic energy. In the usual MP partitioning of the Hamiltonian, where $$F_N$$ is the unperturbed part and $$V_N$$ is the perturbation, the MP2 energy correction is given by

$$\Delta E_{MP2} = \langle \Phi_0 | V_N R^{(0)} V_N | \Phi_0 \rangle$$

where $$R^{(0)} = \sum_{n\neq 0} \dfrac{|\Phi_n\rangle \langle \Phi_n |}{-\Delta_n}$$ is the reduced resolvent operator for the unperturbed Hamiltonian and $$\Delta_n \equiv \Delta_{a_1\ldots a_n}^{i_1\ldots i_n} = \sum_{\gamma = 1}^n (\epsilon_{a_\gamma} - \epsilon_{i_\gamma})$$ is the usual MP energy denominator. Using the diagrammatic rules of many-body perturbation theory, we have the following Goldstone diagrams which I've drawn and evaluated on paper using the usual notation that $$i,j,k,\ldots$$ denote occupied MOs while $$a,b,c,\ldots$$ denote unoccupied MOs (the rules can be found in books like Szabo & Ostlund):

Spin-adaptation diagrammatically is done by associated a weight factor of $$2^l$$ where $$l$$ is the number of closed loops in the Goldstone orbital diagrams. Hence this kind of spin-adaptation can only be done with Goldstones. In the following, we will denote $$h$$ as the number of internally-contracted hole (occupied) lines, $$w$$ as the topological weight of the diagram, and $$s$$ as the sign of the diagram. I'll also show explicitly how this $$2^l$$ factor can be seen by breaking up of spinorbital states $$|p\rangle = |P\rangle \otimes |\sigma_P\rangle$$ where $$P$$ is a spatial orbital index and $$\sigma_P = \pm \frac{1}{2}$$ is the spin projection.

Diagram I: $$w = \frac{1}{2}$$, $$h = 2$$, $$l = 2$$, $$s = (-)^{l+h} = 1$$

$$2^l$$ Rule:

$$\dfrac{1}{2}2^2\sum_{IJAB} \dfrac{\langle IJ|AB\rangle \langle AB|IJ\rangle}{\epsilon_I - \epsilon_A + \epsilon_J - \epsilon_B} = 2\sum_{IJAB} \dfrac{\langle IJ|AB\rangle \langle AB|IJ\rangle}{\epsilon_I - \epsilon_A + \epsilon_J - \epsilon_B}$$

Explicitly integrating out spins:

$$\sum_{ijab} \langle ij|ab\rangle \langle ab|ij\rangle = \sum_{IJAB} \langle IJ|AB\rangle \langle AB|IJ\rangle \sum_{\sigma_A \sigma_B \sigma_I \sigma_J}\langle \sigma_I | \sigma_A \rangle \langle \sigma_J | \sigma_B \rangle \langle \sigma_A | \sigma_I \rangle \langle \sigma_B | \sigma_J \rangle$$

Switching the middle two inner products (since they're scalars)

$$= \sum_{IJAB} \langle IJ|AB\rangle \langle AB|IJ\rangle \sum_{\sigma_A \sigma_B \sigma_I \sigma_J}\langle \sigma_I | \sigma_A \rangle \langle \sigma_A | \sigma_I \rangle \langle \sigma_J | \sigma_B \rangle \langle \sigma_B | \sigma_J \rangle$$

and using resolution of identity in the spin-$$\frac{1}{2}$$ space $$\sum_{\sigma_A} |\sigma_A \rangle \langle \sigma_A | = \sum_{\sigma_B} |\sigma_B \rangle \langle \sigma_B| = 1$$

$$= \sum_{IJAB} \langle IJ|AB\rangle \langle AB|IJ\rangle \sum_{\sigma_I \sigma_J}\langle \sigma_I | \sigma_I \rangle \langle \sigma_J | \sigma_J \rangle = (2)^2 \sum_{IJAB} \langle IJ|AB\rangle \langle AB|IJ\rangle$$

Diagram II: $$w = \frac{1}{2}$$, $$h = 2$$, $$l = 1$$, $$s = (-)^{l+h} = -1$$

$$2^l$$ Rule:

$$-\dfrac{1}{2}2\sum_{IJAB} \dfrac{\langle IJ|BA\rangle \langle AB|IJ\rangle}{\epsilon_I - \epsilon_A + \epsilon_J - \epsilon_B} = -\sum_{IJAB} \dfrac{\langle IJ|BA\rangle \langle AB|IJ\rangle}{\epsilon_I - \epsilon_A + \epsilon_J - \epsilon_B}$$

Explicitly integrating out spins:

$$\sum_{ijab} \langle ij|ba\rangle \langle ab|ij\rangle = \sum_{IJAB} \langle IJ|BA\rangle \langle AB|IJ\rangle \sum_{\sigma_A \sigma_B \sigma_I \sigma_J} \langle \sigma_I | \sigma_B \rangle \langle \sigma_J | \sigma_A \rangle \langle \sigma_A | \sigma_I \rangle \langle \sigma_B | \sigma_J \rangle$$

and using resolution of identity again to set $$\sum_{\sigma_A} |\sigma_A \rangle \langle \sigma_A | = 1$$

$$= \sum_{IJAB} \langle IJ|BA\rangle \langle AB|IJ\rangle \sum_{\sigma_B \sigma_I \sigma_J} \langle \sigma_I | \sigma_B \rangle \langle \sigma_J | \sigma_I \rangle \langle \sigma_B | \sigma_J \rangle$$

using $$\langle \sigma_I | \sigma_J \rangle = \delta_{IJ}$$

$$= \sum_{IJAB} \langle IJ|BA\rangle \langle AB|IJ\rangle \sum_{\sigma_B \sigma_I} \langle \sigma_I | \sigma_B \rangle \langle \sigma_B | \sigma_I \rangle$$

using $$\sum_{\sigma_B} |\sigma_B \rangle \langle \sigma_B | = 1$$

$$= \sum_{IJAB} \langle IJ|BA\rangle \langle AB|IJ\rangle \sum_{\sigma_I} \langle \sigma_I | \sigma_I \rangle = 2\sum_{IJAB} \langle IJ|BA\rangle$$

So we can see that the $$2^l$$ rule really works (and of course, there is a more general proof that covers all cases)

Thus, the total MP2 correction in closed-shell singlet spin-adapted form is: $$\Delta E_{MP2} = \sum_{IJAB} \dfrac{2\langle IJ| AB\rangle \langle AB | IJ\rangle - \langle IJ | BA \rangle \langle AB | IJ \rangle }{\epsilon_I - \epsilon_A + \epsilon_J - \epsilon_B}$$

Only in the case of real orbitals (usual case), we can write this as the more compact $$\Delta E_{MP2} = \sum_{IJAB} \dfrac{2\langle IJ| AB\rangle ^2 - \langle IJ | BA \rangle \langle AB | IJ \rangle }{\epsilon_I - \epsilon_A + \epsilon_J - \epsilon_B}$$