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

Question

$\newcommand{\Ket}[1]{\left|#1\right>}$ $\newcommand{\Bra}[1]{\left<#1\right|}$ $\newcommand{\BraKet}[2] { {\left<#1} \left|#2 \right>}$
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.

Answer 1

$\newcommand{\Ket}[1]{\left|#1\right>}$ $\newcommand{\Bra}[1]{\left<#1\right|}$ $\newcommand{\BraKet}[2] { {\left<#1} \left|#2 \right>}$
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.

Answer 2

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