Skip to main content

Toggle Sidebar Navigation Menu

Quantum Chemistry Method Development and Simulations

• Home
• Research
• People
• Publications
• Software
• Tutorials
• Funding
• Photos
• Group Area

Powered by Jekyll & Tao

Many-Body Perturbation Theory

Jupyter Notebook Tutorials: Quantum Chemistry

---

1. Hartree Fock Theory
2. Second Quantization
3. Many-Body Perturbation Theory
4. Configuration Interaction with Single Excitations
5. Static Response Properties

---

Many-Body Perturbation Theory

The Jupyter notebook for this tutorial can be found here.

Many-body perturbation theory (MBPT) is a non-variational approach based on non-degenerate perturbation theory that can be used to find approximate solutions to the electronic Schrödinger equation. This project provides an overview of the theoretical aspects of MBPT and a derivation of energy and wave function expressions for the first few orders, aided by the p$^\dagger$q package for manipulating strings of second-quantized operators. For a review of key concepts in non-degenerate perturbation theory, refer to these notes. For an introduction to the p$^\dagger$q package, refer to this tutorial.

The Hamiltonian

We begin by partitioning the Hamiltonian for a many-electron system into a "simple" Hamiltonian for which solutions can be determined and a "small" perturbation to that simple Hamiltonian. Here, we use the partitioning

$$\begin{align} \hat{H} = \hat{f} + \lambda \hat{v} \end{align}$$

where the zeroth-order Hamiltonian, $\hat{f}$, is the Fock operator, the perturbation, $\hat{v}$, is called the fluctuation potential, and $\lambda$ is the perturbation strength. Recall from the Hartree-Fock tutorial that the second-quantized Fock operator is

$$\begin{align} \hat{f} = \sum_{\mu\nu} f_{\mu\nu}\hat{a}^\dagger_\mu \hat{a}_\nu \end{align}$$

with

$$\begin{align} f_{\mu\nu} = h_{\mu\nu} + J_{\mu\nu} - K_{\mu\nu} \end{align}$$

Here, $\mu$ and $\nu$ are labels for the atomic orbital (AO) basis functions, and $h_{\mu\nu}$, $J_{\mu\nu}$, and $K_{\mu\nu}$ are elements of the core Hamiltonian, Coulomb, and exchange matrices, respectively. In the canonical molecular orbital (MO) basis determined via the Hartree-Fock-Roothaan procedure, we have

$$\begin{align} \hat{f} = \sum_{pq} f_{pq}\hat{a}^\dagger_p \hat{a}_q \end{align}$$

with

$$\begin{align} f_{pq} = h_{pq} + \sum_i \langle pi||qi \rangle = \epsilon_p \delta_{pq} \end{align}$$

Here, $p$ and $q$ are general MO labels spanning the occupied and virtual orbital spaces, and $i$ is an occupied MO label. The symbol $h_{pq}$ is again an element of the core Hamiltonian matrix, $\langle pi||qi \rangle$ is an antisymmetrized two-electron repulsion integral in physicists' notation, and $\epsilon_p$ is the energy of MO $p$. Note that the Fock matrix is diagonal in the MO basis, so the Fock operator is simply

$$\begin{align} \hat{f} = \sum_p \epsilon_p \hat{a}^\dagger_p \hat{a}_p \end{align}$$

The full second-quantized Hamiltonian represented in the MO basis is

$$\begin{align} \hat{H} = \sum_{pq} h_{pq} + \frac{1}{4} \sum_{pqrs} \hat{a}^\dagger_p \hat{a}_q^\dagger \hat{a}_s\hat{a}_r \langle pq||rs\rangle \end{align}$$

Given the definition of the Fock operator above, the fluctuation potential must then be

$$\begin{align} \hat{v} = \frac{1}{4} \sum_{pqrs} \hat{a}^\dagger_p \hat{a}_q^\dagger \hat{a}_s\hat{a}_r \langle pq||rs\rangle - \sum_{pq} \hat{a}^\dagger_p \hat{a}_q \sum_i \langle pi||qi \rangle \end{align}$$

The Perturbation Expansion of the Wave Function and Energy

Recall that, in non-degenerate perturbation theory, the wave function and energy are expanded in powers of the perturbation strength as

$$\begin{align} |\psi_n\rangle = |\psi_n^{(0)}\rangle + \lambda |\psi_n^{(1)}\rangle + \lambda ^2|\psi_n^{(2)}\rangle + \lambda ^3|\psi_n^{(3)}\rangle + ... \end{align}$$

and

$$\begin{align} E_n = E_n^{(0)} + \lambda E_n^{(1)} + \lambda ^2 E_n^{(2)} + \lambda^3 E_n^{(3)} + ... \end{align}$$

where $|\psi_n^{(n)}\rangle$ is the $p$-th order correction to the wave function and $E_n^{(p)}$ is the $p$-th order correction to the energy. In MBPT, these quantities are determined order by order, beginning with the zeroth-order ($p = 0$) solution to the problem.

What are the zeroth-order wave functions $|\psi_n^{(0)}\rangle$? Is the Hartree-Fock wave function ($|\Phi_0\rangle$) a suitable zeroth-order wave function? The eigenfunctions and eigenvalues of the Fock operator are the MOs and corresponding MO energies that satisfy

$$\begin{align} \hat{f} \phi_p = \epsilon_p \phi_p \end{align}$$

The Hatree-Fock wave function is a determinant of these orbitals, i.e.,

$$ \Phi_0 ({\bf x}_1, {\bf x}_2, ..., {\bf x}_N) = \frac{1}{\sqrt{N!}}\begin{vmatrix} \phi_1({\bf x}_1) & \phi_2({\bf x}_1) & ... & \phi_N({\bf x}_1) \\ \phi_1({\bf x}_2) & \phi_2({\bf x}_2) & ... & \phi_N({\bf x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_1({\bf x}_N) & ... & ... & \phi_N({\bf x}_N) \end{vmatrix} $$

where $N$ is the number of electrons in the system, and ${\bf x}_i$ refers to the spin and spatial coordinates of electron $i$. An antisymmetrized product of MOs such as this Slater determinant is also an eigenfunction of $\hat{f}$ that satisfies

$$\begin{align} \hat{f}|\Phi_0\rangle = E^{(0)}_0 |\Phi_0\rangle \end{align}$$

with

$$\begin{align} E^{(0)}_0 = \sum_i \epsilon_i \end{align}$$

where the sum is restricted to include only those orbitals that are occupied in $|\Phi_0\rangle$. Throughout this tutorial, the labels $i$, $j$, $k$, and $l$ refer to MOs that are occupied in the Hartree-Fock configuration, while $a$, $b$, $c$, and $d$ refer to virtual MOs.

In MBPT, we choose the zeroth-order wave functions to be all possible $N$-electron Slater determinants, including the Hartree-Fock configuration, i.e.,

$$\begin{align} |\psi_0^{(0)}\rangle &= |\Phi_0\rangle \\ |\psi_{a,i}^{(0)}\rangle &= |\Phi_i^a \rangle = \hat{a}^\dagger_a \hat{a}_i |\Phi_0\rangle \\ |\psi_{ab,ij}^{(0)}\rangle &= |\Phi_{ij}^{ab} \rangle = \hat{a}^\dagger_a \hat{a}^\dagger_b \hat{a}_j \hat{a}_i |\Phi_0\rangle \\ \text{etc.} \nonumber \end{align}$$

Here, we have introduced the notation $|\Phi_i^a\rangle$, $|\Phi_{ij}^{ab}\rangle$, etc. to refer to determinants of orbitals that differ from the Hartree-Fock determinant by a single electron transition, a double electron transition, etc. In a basis of orthonormal MOs, the wave functions for these excited configurations are orthogonal the Hartree-Fock wave function, meaning that they satisfy

$$\begin{align} \langle \Phi_i^a | \Phi_0 \rangle &= 0 \\ \langle \Phi_{ij}^{ab} | \Phi_0 \rangle &= 0 \\ \text{etc.} \nonumber \end{align}$$

which is consistent with the usual requirement in non-degenerate perturbation theory that the zeroth-order wave functions form an orthonormal set.

The zeroth-order energies corresponding to these zeroth-order wave functions are sums of the energies of the orbitals that are occupied within each configuration. For example, the zeroth-order energy for a singly-excited configuration, $|\Phi_i^a\rangle$ is

$$\begin{align} E^{(0)}_{a,i} = \sum_{k\neq i} \epsilon_k + \epsilon_a \end{align}$$

Let us confirm our understanding of the zeroth-order energies using the p$^\dagger$q package. The following code will evaluate the expectation value of the Fock operator with respect to the Hartree-Fock determinant (the Fermi vacuum).

import pdaggerq

pq = pdaggerq.pq_helper("fermi")

pq.add_operator_product(1.0, ['f'])
pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

['+1.00', 'f(i,i)']

Great! This is exactly what we were expecting. The expectation value of the Fock operator is equal to the sum over the diagonal elements of the Fock matrix (the orbital energies), and the sum is restricted to the occupied orbital space.

For the singly excited functions, we can do something similar, but we need to be aware of a key internal assumption made by the p$^\dagger$q package, that repeated labels imply a summation. As such, an expression like

$$\begin{align} E^{(0)}_{a,i} = \langle \Phi^a_i | \hat{f} | \Phi_i^a \rangle = \langle \Phi_0 | \hat{a}^\dagger_i \hat{a}_a \hat{f} \hat{a}^\dagger_a \hat{a}_i | \Phi_{0}\rangle \end{align}$$

would be treated as a contraction over labels $i$ and $a$. We can avoid this technical issue by evaluating a matrix element involving different labels in the bra and ket states as

import pdaggerq

pq = pdaggerq.pq_helper("fermi")

pq.set_left_operators([['a*(i)', 'a(a)']])
pq.set_right_operators([['a*(b)', 'a(j)']])

pq.add_operator_product(1.0, ['f'])
pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

['+1.00', 'd(a,b)', 'd(i,j)', 'f(k,k)']
['-1.00', 'd(a,b)', 'f(j,i)']
['+1.00', 'd(i,j)', 'f(a,b)']

Here, the notation 'd(a,b)' refers to a Kronecker delta function, $\delta_{ab}$. Now, we can think about the case where $i = j$ and $a = b$ and see that the result is

$$\begin{align} \langle \Phi_0 | \hat{a}^\dagger_i \hat{a}_a \hat{f} \hat{a}^\dagger_a \hat{a}_i | \Phi_{0}\rangle &= \sum_k f_{kk} - f_{ii} + f_{aa} \nonumber \\ &= \sum_k \epsilon_{k} - \epsilon_{i} + \epsilon_{a} \nonumber \\ &= \sum_{k\neq i} \epsilon_{k} + \epsilon_{a} \nonumber \\ &= E^{(0)}_{a,i} \end{align}$$

which is exactly what we were expecting.

The First-Order Correction to the Energy

Now, we proceed by determining the corrections to the energy and wave function, using the Hartree-Fock wave function as the target zeroth-order state. Let us consider the first-order part of the Schrödinger equation

$$\begin{align} \hat{f} |\psi_n^{(1)}\rangle + \hat{v} | \psi_n^{(0)} \rangle = E_n^{(0)} |\psi_n^{(1)}\rangle + E_n^{(1)} |\psi_n^{(0)}\rangle \end{align}$$

We take $n = 0$ and $| \psi_0^{(0)} \rangle = |\Phi_0\rangle$, which then gives us

$$\begin{align} \hat{f} |\psi_0^{(1)}\rangle + \hat{v} | \Phi_0 \rangle = E_0^{(0)} |\psi_0^{(1)}\rangle + E_0^{(1)} |\Phi_0 \rangle \end{align}$$

We project this equation onto $\langle \Phi_0|$ to obtain

$$\begin{align} \langle \Phi_0 | \hat{f} |\psi_0^{(1)}\rangle + \langle \Phi_0 | \hat{v} | \Phi_0 \rangle = E_0^{(0)} \langle \Phi_0 |\psi_0^{(1)}\rangle + E_0^{(1)} \langle \Phi_0 |\Phi_0 \rangle \\ \end{align}$$

which simplifies to

$$\begin{align} \langle \Phi_0 | \hat{v} | \Phi_0 \rangle = E_0^{(1)} \end{align}$$

because

$$\begin{align} \langle \Phi_0 | \hat{f} |\psi_0^{(1)}\rangle = E_0^{(0)} \langle \Phi_0 | \psi_0^{(1)}\rangle = 0 \end{align}$$ and $$\begin{align} \langle \Phi_0 | \Phi_0 \rangle = 1 \end{align}$$

Now, let us use p$^\dagger$q to evaluate the first-order correction to the energy. The following code evaluates the expectation value of the fluctuation potential with respect to the Fermi vacuum, which is the target zeroth-order state.

pq = pdaggerq.pq_helper("fermi")

pq.add_operator_product(1.0, ['v'])
pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

['-0.50', '<j,i||j,i>']

We can see that, to first-order, the energy for state $|\Phi_0\rangle$ is

$$\begin{align} E_0 &= E_0^{(0)} + E_1^{(0)} \nonumber \\ &= \sum_i \epsilon_i - \frac{1}{2}\sum_{ij} \langle ij||ij \rangle \end{align}$$

which is simply the Hartree-Fock energy! Remarkable!

The First-Order Correction to the Wave Function

The first-order correction to the wave function, $|\psi_0^{(1)}\rangle$, can be expanded in terms of zeroth-order wave functions as

$$\begin{align} |\psi_0^{(1)}\rangle = \sum_{k \neq 0} c^{(1)}_k |\psi_k^{(0)}\rangle \end{align}$$

where the restriction on the sum ensures that the first-order correction to the wave function for state 0 is orthogonal to $|\psi_0^{(0)}\rangle$. If we recall that our zeroth-order wave functions are just Slater determinants, we can recast the equation above as

$$\begin{align} |\psi_0^{(1)}\rangle = \sum_{i,a} c_{a,i}^{(1)} |\Phi_i^a\rangle + \frac{1}{4} \sum_{ij,ab} c_{ab, ij}^{(1)} |\Phi_{ij}^{ab}\rangle + \frac{1}{36} \sum_{ijk, abc} c_{abc, ijk}^{(1)} |\Phi_{ijk}^{abc}\rangle + ... \end{align}$$

The first-order expansion coefficients in this expression are defined by the projections

$$\begin{align} c_{a,i}^{(1)} &= \langle \Phi_i^a|\psi_0^{(1)}\rangle \\ c_{ab,ij}^{(1)} &= \langle \Phi_{ij}^{ab}|\psi_0^{(1)}\rangle \\ \text{etc.}\nonumber \end{align}$$

Now, consider again the first-order part of the Schrödinger equation

$$\begin{align} \hat{f} |\psi_0^{(1)}\rangle + \hat{v} | \Phi_0 \rangle = E_0^{(0)} |\psi_0^{(1)}\rangle + E_0^{(1)} |\Phi_0 \rangle \end{align}$$

We project this equation onto a zeroth-order wave function that is not equal to $\langle\Phi_0|$, e.g., $\langle \Phi_{ij...}^{ab...}|$, to obtain

$$\begin{align} \langle \Phi_{ij...}^{ab...} | \hat{f} |\psi_0^{(1)}\rangle + \langle \Phi_{ij...}^{ab...} | \hat{v} | \Phi_0 \rangle = E_0^{(0)} \langle \Phi_{ij...}^{ab...} |\psi_0^{(1)}\rangle + E_0^{(1)} \langle \Phi_{ij...}^{ab...} |\Phi_0 \rangle \\ \end{align}$$

which simplifies to

$$\begin{align} E_{ab...,ij...}^{(0)} c_{ab...,ij...}^{(1)} + \langle \Phi_{ij...}^{ab...} | \hat{v} | \Phi_0 \rangle = E_0^{(0)} c_{ab...,ij...}^{(1)} \end{align}$$

because

$$\begin{align} \langle \Phi_{ij...}^{ab...} | \hat{f} |\psi_0^{(1)}\rangle = E_{ab...,ij...}^{(0)} \langle \Phi_{ij...}^{ab...} |\psi_0^{(1)} \rangle = E_{ab...,ij...}^{(0)} c_{ab...,ij...}^{(1)} \end{align}$$ and $$\begin{align} \langle \Phi_{ij...}^{ab...} |\Phi_0 \rangle = 0 \end{align}$$

Solving for the first-order amplitudes, we have

$$\begin{align} c_{a,i}^{(1)} &= \frac{\langle \Phi_i^a | \hat{v} | \Phi_0 \rangle}{E_0^{(0)} - E_{a,i}^{(0)}} \\ c_{ab,ij}^{(1)} &= \frac{\langle \Phi_{ij}^{ab} | \hat{v} | \Phi_0 \rangle}{E_0^{(0)} - E_{ab,ij}^{(0)}} \\ \text{etc.}\nonumber \end{align}$$

or

$$\begin{align} c_{a,i}^{(1)} &= \frac{\langle \Phi_i^a | \hat{v} | \Phi_0 \rangle}{\epsilon_i - \epsilon_a} \\ c_{ab,ij}^{(1)} &= \frac{\langle \Phi_{ij}^{ab} | \hat{v} | \Phi_0 \rangle}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b} \\ \text{etc.}\nonumber \end{align}$$

We can use p$^\dagger$q to evaluate the numerators in these expressions using the following code.

pq = pdaggerq.pq_helper("fermi")

print('[E_0 - E(a,i)]c(a,i) = <a,i|v|0>:')
print('')

pq.set_left_operators([['a*(i)', 'a(a)']])
pq.add_operator_product(1.0, ['v'])
pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

print('')
print('[E_0 - E(ab,ij)]c(ab,ij) = <ab,ij|v|0>:')
print('')

pq.set_left_operators([['a*(i)', 'a*(j)', 'a(b)', 'a(a)']])
pq.add_operator_product(1.0, ['v'])
pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

print('')
print('[E_0 - E(abc,ijk)]c(abc,ijk) = <abc,ijk|v|0>:')
print('')

pq.set_left_operators([['a*(i)', 'a*(j)', 'a*(k)', 'a(c)', 'a(b)', 'a(a)']])
pq.add_operator_product(1.0, ['v'])
pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

[E_0 - E(a,i)]c(a,i) = <a,i|v|0>:


[E_0 - E(ab,ij)]c(ab,ij) = <ab,ij|v|0>:

['+1.00', '<a,b||i,j>']

[E_0 - E(abc,ijk)]c(abc,ijk) = <abc,ijk|v|0>:

From these expressions, it appears that

$$\begin{align} c^{(1)}_{a,i} &= 0 \\ c^{(1)}_{ab,ij} &= \frac{\langle ab || ij \rangle }{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b} \\ c^{(1)}_{abc,ijk} &= 0 \\ \text{etc.} \nonumber \end{align}$$

The fact that the first-order singles amplitudes are zero is a consequence of Brillouin's theorem. We should have anticipated that the first-order triples amplitudes (and first-order quadruples amplitudes, etc.) would be zero because the fluctuation potential is a two-body operator, and the Slater-Condon rules tell us that matrix elements of two body operators involving determinants that differ by more than two electrons will always be zero.

Our final expression for the first-order correction to the wave function is

$$\begin{align} |\psi_0^{(1)}\rangle = \frac{1}{4} \sum_{ij,ab} |\Phi_{ij}^{ab}\rangle \frac{\langle ab || ij \rangle }{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b} \end{align}$$

The Second-Order Correction to the Energy

We found earlier that the energy to first-order in MBPT is simply the Hartree-Fock energy. As such, the first meaningful correction to the energy is obtained at second order. To find, $E_0^{(2)}$, let us consider the second-order part of the Schrödinger equation

$$\begin{align} \hat{f} |\psi_0^{(2)}\rangle + \hat{v} |\psi_0^{(1)} \rangle = E_0^{(0)} |\psi_0^{(2)}\rangle + E_0^{(1)}|\psi_0^{(1)}\rangle + E_0^{(2)} |\Phi_0 \rangle \end{align}$$

When seeking the first-order correction to the energy, we projected the first-order Schrödinger equation onto $\langle \Phi_0|$. Doing the same for the second-order Schrödinger equation leads to

$$\begin{align} \langle \Phi_0 | \hat{f} |\psi_0^{(2)}\rangle + \langle \Phi_0 | \hat{v} |\psi_0^{(1)} \rangle = E_0^{(0)} \langle \Phi_0|\psi_0^{(2)}\rangle + E_0^{(1)}\langle \Phi_0|\psi_0^{(1)}\rangle + E_0^{(2)} \langle \Phi_0|\Phi_0 \rangle \end{align}$$

which simplifies to

$$\begin{align} \langle \Phi_0 | \hat{v} | \psi^{(1)}_0 \rangle = E_0^{(2)} \end{align}$$

because

$$\begin{align} \langle \Phi_0 | \hat{f} |\psi_0^{(2)}\rangle = E_0^{(0)} \langle \Phi_0 | \psi_0^{(2)}\rangle = 0 \end{align}$$ and, of course, $$\begin{align} \langle \Phi_0 | \Phi_0 \rangle = 1 \end{align}$$

Inserting the definition of the first-order correction to the wave function gives us

$$\begin{align} E_0^{(2)} = \frac{1}{4} \sum_{ij,ab} \langle \Phi_0 | \hat{v}|\Phi_{ij}^{ab}\rangle \frac{\langle ab || ij \rangle }{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b} \end{align}$$

We have already evaluated the matrix element, $\langle \Phi_0 | \hat{v}|\Phi_{ij}^{ab}\rangle$, so the second-order correction to the energy is actually given by

$$\begin{align} E_0^{(2)} = \frac{1}{4} \sum_{ij,ab} \frac{|\langle ab || ij \rangle |^2}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b} \end{align}$$

Note that, in this expression, the labels $i$, $j$, $a$, and $b$ refer to spin orbitals. If we are working in a spatial orbital basis, we must explicitly account for the spin of each orbital label. In this case, we would have

$$\begin{align} E_0^{(2)} = \frac{1}{4} \sum_\sigma \sum_{ij,ab} \frac{|\langle a_\sigma b_\sigma || i_\sigma j_\sigma \rangle |^2}{\epsilon_{i_\sigma} + \epsilon_{j_\sigma} - \epsilon_{a_\sigma} - \epsilon_{b_\sigma}} + \frac{1}{2} \sum_{\sigma\neq\tau} \sum_{ij,ab} \frac{|\langle a_\sigma b_\tau | i_\sigma j_\tau \rangle |^2}{\epsilon_{i_\sigma} + \epsilon_{j_\tau} - \epsilon_{a_\sigma} - \epsilon_{b_\tau}} \end{align}$$

where $\sigma, \tau \in \{\alpha, \beta\}$. Note the factor of $\frac{1}{2}$ and the lack of antisymmetrization on the opposite-spin term, which arise from spin-symmetry considerations, i.e., $\langle a_\sigma b_\tau || i_\sigma j_\tau \rangle = \langle a_\sigma b_\tau | i_\sigma j_\tau \rangle$ and $\langle a_\sigma b_\tau || i_\tau j_\sigma \rangle = -\langle a_\sigma b_\tau | j_\sigma i_\tau \rangle$.

Numerical Procedure for Evaluating $E_0^{(2)}$

To evaluate the second-order energy numerically, we can use the Psi4 package (or any other package that gives us access to the required integrals). The following code runs a Hartree-Fock calculation on H$_2$O$^+$ and extracts the $\alpha$- and $\beta$-spin MO coefficient matrices and energies. Note that some of the keyword options are chosen to tighten convergence and to avoid the use of the density fitting approximation. Eventually, we will compare our second-order energy (and, later on, the third-order energy) to the same quantity evaluated using Psi4, and we want to avoid any sources of possible disagreement.

import psi4
import numpy as np

mol = psi4.geometry("""
1 2
O
H 1 1.0
H 1 1.0 2 104.5
symmetry c1
""")

options = {'basis' : 'sto-3g', 
           'scf_type': 'pk', 
           'mp2_type' : 'conv', 
           'mp_type' : 'conv',
           'e_convergence' : 1e-10,
           'd_convergence' : 1e-10,
           'reference' : 'uhf'}

psi4.set_options(options)
psi4.set_output_file('output.dat') 

e_hf, wfn = psi4.energy('scf', return_wfn = True)

# get MO coefficients
Ca = wfn.Ca()
Cb = wfn.Cb()

# orbital energies
epsilon_a = wfn.epsilon_a().to_array()
epsilon_b = wfn.epsilon_b().to_array()

We can now use Psi4's MintsHelper object to evaluate the electron repulsion integrals in the MO basis. We will use the spatial orbital expression for $E_0^{(2)}$, which means that we need to consider multiple spin blocks of the electron repulsion integral (ERI) tensor.

# molecular integrals helper
mints = psi4.core.MintsHelper(wfn.basisset())

# ERIs in chemists' notation
g_aaaa = np.array(mints.mo_eri(Ca, Ca, Ca, Ca))
g_bbbb = np.array(mints.mo_eri(Cb, Cb, Cb, Cb))
g_aabb = np.array(mints.mo_eri(Ca, Ca, Cb, Cb))

# ERIs in physicists' notation <pq|rs> = (pr|qs)
g_aaaa = g_aaaa.transpose(0, 2, 1, 3)
g_bbbb = g_bbbb.transpose(0, 2, 1, 3)
g_abab = g_aabb.transpose(0, 2, 1, 3)

# antisymmetrized ERIs in physicists' notation <pq||rs> = <pq|rs> - <pq|sr>
g_aaaa -= g_aaaa.transpose(0, 1, 3, 2)
g_bbbb -= g_bbbb.transpose(0, 1, 3, 2)

The following code evaluates $E_0^{(2)}$ and compares the result to the correlation energy from an MP2 calculation carried out using Psi4. Second-order MBPT is commonly referred to as second-order Møller-Plesset perturbation theory, or MP2.

# number of electrons, basis functions
nmo = wfn.nmo()
nalpha = wfn.nalpha()
nbeta = wfn.nbeta()

e2 = 0.0
# aaaa
for i in range (0, nalpha):
    for j in range (0, nalpha):
        for a in range (nalpha, nmo):
            for b in range (nalpha, nmo):
                denom = epsilon_a[i] + epsilon_a[j] - epsilon_a[a] - epsilon_a[b]
                e2 += 0.25 * g_aaaa[i, j, a, b]**2 / denom

# bbbb
for i in range (0, nbeta):
    for j in range (0, nbeta):
        for a in range (nbeta, nmo):
            for b in range (nbeta, nmo):
                denom = epsilon_b[i] + epsilon_b[j] - epsilon_b[a] - epsilon_b[b]
                e2 += 0.25 * g_bbbb[i, j, a, b]**2 / denom


# abab and baba (note the factor of two difference from the expression above)
for i in range (0, nalpha):
    for j in range (0, nbeta):
        for a in range (nalpha, nmo):
            for b in range (nbeta, nmo):
                denom = epsilon_a[i] + epsilon_b[j] - epsilon_a[a] - epsilon_b[b]
                e2 += g_abab[i, j, a, b]**2 / denom
                
e_mp2 = psi4.energy('mp2')
e2_psi4 = e_mp2 - e_hf

print('')
print('my E(2) =   %20.12f Eh' % (e2))
print('psi4 E(2) = %20.12f Eh' % (e2_psi4))
print('')

my E(2) =        -0.029933352948 Eh
psi4 E(2) =      -0.029933352948 Eh

We can see that our second-order energy is numerically equivalent to that from Psi4. Fantastic!

The Third-Order Correction to the Energy

To find, $E_0^{(3)}$, let us consider the third-order part of the Schrödinger equation

$$\begin{align} \hat{f} |\psi_0^{(3)}\rangle + \hat{v} |\psi_0^{(2)} \rangle = E_0^{(0)} |\psi_0^{(3)}\rangle + E_0^{(1)}|\psi_0^{(2)}\rangle + E_0^{(2)}|\psi_0^{(1)}\rangle + E_0^{(3)} |\Phi_0 \rangle \end{align}$$

We project this equation onto $\langle \Phi_0|$ to obtain

$$\begin{align} \langle \Phi_0 | \hat{f} |\psi_0^{(3)}\rangle + \langle \Phi_0 | \hat{v} |\psi_0^{(2)} \rangle = E_0^{(0)} \langle \Phi_0 |\psi_0^{(3)}\rangle + E_0^{(1)}\langle \Phi_0 |\psi_0^{(2)}\rangle + E_0^{(2)}\langle \Phi_0 |\psi_0^{(1)}\rangle + E_0^{(3)} \langle \Phi_0 |\Phi_0 \rangle\end{align}$$

which simplifies to

$$\begin{align} \langle \Phi_0 | \hat{v} | \psi^{(2)}_0 \rangle = E_0^{(3)} \end{align}$$

because

$$\begin{align} \langle \Phi_0 | \psi_0^{(2)}\rangle &= 0 \\ \langle \Phi_0 | \hat{f} |\psi_0^{(3)}\rangle &= E_0^{(0)} \langle \Phi_0 | \psi_0^{(3)}\rangle = 0 \end{align}$$ and, yet again, $$\begin{align} \langle \Phi_0 | \Phi_0 \rangle = 1 \end{align}$$

Oh, no! To evaluate the third-order energy, we need the second-order correction to the wave function! Fortunately, the situation is not as terrible as it seems because we only need the doubles part of the second-order correction to the wave function. The reason is that the matrix element $\langle \Phi_0 | \hat{v} | \psi^{(2)}_0 \rangle$ is only non-zero when the bra states correspond to doubly-excited determinants, as we saw earlier when evaluating the first-order correction to the wave function. This result does not imply that the second-order correction to the wave function contains only doubles terms; that is not the case. Rather, all we have shown is that the third-order correction to the energy can be determined from this subset of the second-order correction to the wave function, i.e.,

$$\begin{align} E_0^{(3)} = \frac{1}{4} \sum_{ij,ab} \langle \Phi_0 | \hat{v} | \Phi_{ij}^{ab}\rangle c^{(2)}_{ab,ij} = \frac{1}{4} \sum_{ij,ab} \langle ij||ab \rangle c^{(2)}_{ab,ij} \end{align}$$

To obtain the second-order doubles amplitudes, we project the second-order part of the Schrödinger equation onto $\langle \Phi_{ij}^{ab}|$ giving

$$\begin{align} \langle \Phi_{ij}^{ab} | \hat{f} |\psi_0^{(2)}\rangle + \langle \Phi_{ij}^{ab} | \hat{v} |\psi_0^{(1)} \rangle = E_0^{(0)} \langle \Phi_{ij}^{ab}|\psi_0^{(2)}\rangle + E_0^{(1)}\langle \Phi_{ij}^{ab}|\psi_0^{(1)}\rangle + E_0^{(2)} \langle \Phi_{ij}^{ab}|\Phi_0 \rangle \end{align}$$

or

$$\begin{align} E_{ab,ij}^{(0)} c^{(2)}_{ab,ij} + \langle \Phi_{ij}^{ab} | \hat{v} |\psi_0^{(1)} \rangle = E_0^{(0)}c^{(2)}_{ab,ij} + E_0^{(1)}c^{(1)}_{ab,ij} \end{align}$$

Solving for the second-order doubles amplitudes, we have

$$\begin{align} c^{(2)}_{ab,ij} = \frac{\langle \Phi_{ij}^{ab} | \hat{v} |\psi_0^{(1)} \rangle - E_0^{(1)}c^{(1)}_{ab,ij}}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b} \end{align}$$

Since we have alread determined $E_0^{(1)}$ and $c^{(1)}_{ab,ij}$, the only remaining challenge is to evaluate

$$\begin{align} \langle \Phi_{ij}^{ab} | \hat{v} |\psi_0^{(1)} \rangle = \frac{1}{4} \sum_{kl,cd} \langle \Phi_{ij}^{ab} | \hat{v} | \Phi_{kl}^{cd} \rangle c^{(1)}_{cd,kl} \end{align}$$

The following code evaluates this quantity using p$^\dagger$q. To simplify things, we will make use of the built-in operator 'r2', the action of which on the Hartree-Fock wave function can represent the first-order wave function. We will introduce spin labels, which can be done by passing a dictionary of spins to the pq.strings() function.

pq = pdaggerq.pq_helper('fermi')

pq.set_left_operators([['a*(i)', 'a*(j)', 'a(b)', 'a(a)']])
pq.set_right_operators([['r2']])
pq.add_operator_product(1.0, ['v'])
pq.simplify()

print('')
print('# aa terms')
print('')

spin_labels = {
    'a' : 'a',
    'b' : 'a',
    'i' : 'a',
    'j' : 'a'
}

terms = pq.strings(spin_labels = spin_labels)
for term in terms:
    print(term)

print('')
print('# bb terms')
print('')

spin_labels = {
    'a' : 'b',
    'b' : 'b',
    'i' : 'b',
    'j' : 'b'
}

terms = pq.strings(spin_labels = spin_labels)
for term in terms:
    print(term)

print('')
print('# ab terms')
print('')

spin_labels = {
    'a' : 'a',
    'b' : 'b',
    'i' : 'a',
    'j' : 'b'
}

terms = pq.strings(spin_labels = spin_labels)
for term in terms:
    print(term)

# aa terms

['-0.50', '<l,k||l,k>_aaaa', 'r2_aaaa(a,b,i,j)']
['-0.50', '<l,k||l,k>_abab', 'r2_aaaa(a,b,i,j)']
['-0.50', '<k,l||k,l>_abab', 'r2_aaaa(a,b,i,j)']
['-0.50', '<l,k||l,k>_bbbb', 'r2_aaaa(a,b,i,j)']
['+0.50', '<l,k||i,j>_aaaa', 'r2_aaaa(a,b,l,k)']
['+1.00', 'P(i,j)', 'P(a,b)', '<k,a||c,j>_aaaa', 'r2_aaaa(c,b,i,k)']
['-1.00', 'P(i,j)', 'P(a,b)', '<a,k||j,c>_abab', 'r2_abab(b,c,i,k)']
['+0.50', '<a,b||c,d>_aaaa', 'r2_aaaa(c,d,i,j)']

# bb terms

['-0.50', '<l,k||l,k>_aaaa', 'r2_bbbb(a,b,i,j)']
['-0.50', '<l,k||l,k>_abab', 'r2_bbbb(a,b,i,j)']
['-0.50', '<k,l||k,l>_abab', 'r2_bbbb(a,b,i,j)']
['-0.50', '<l,k||l,k>_bbbb', 'r2_bbbb(a,b,i,j)']
['+0.50', '<l,k||i,j>_bbbb', 'r2_bbbb(a,b,l,k)']
['-1.00', 'P(i,j)', 'P(a,b)', '<k,a||c,j>_abab', 'r2_abab(c,b,k,i)']
['+1.00', 'P(i,j)', 'P(a,b)', '<k,a||c,j>_bbbb', 'r2_bbbb(c,b,i,k)']
['+0.50', '<a,b||c,d>_bbbb', 'r2_bbbb(c,d,i,j)']

# ab terms

['-0.50', '<l,k||l,k>_aaaa', 'r2_abab(a,b,i,j)']
['-0.50', '<l,k||l,k>_abab', 'r2_abab(a,b,i,j)']
['-0.50', '<k,l||k,l>_abab', 'r2_abab(a,b,i,j)']
['-0.50', '<l,k||l,k>_bbbb', 'r2_abab(a,b,i,j)']
['+0.50', '<l,k||i,j>_abab', 'r2_abab(a,b,l,k)']
['+0.50', '<k,l||i,j>_abab', 'r2_abab(a,b,k,l)']
['-1.00', '<a,k||c,j>_abab', 'r2_abab(c,b,i,k)']
['-1.00', '<k,b||c,j>_abab', 'r2_aaaa(c,a,i,k)']
['+1.00', '<k,b||c,j>_bbbb', 'r2_abab(a,c,i,k)']
['+1.00', '<k,a||c,i>_aaaa', 'r2_abab(c,b,k,j)']
['-1.00', '<a,k||i,c>_abab', 'r2_bbbb(c,b,j,k)']
['-1.00', '<k,b||i,c>_abab', 'r2_abab(a,c,k,j)']
['+0.50', '<a,b||c,d>_abab', 'r2_abab(c,d,i,j)']
['+0.50', '<a,b||d,c>_abab', 'r2_abab(d,c,i,j)']

Here, 'P(i,j)' represents a permutation operator, and we interpret 'r2(a,b,i,j)' as a first-order doubles amplitude. We can see that these expressions are a bit more complicated than the ones we generated earlier, so it might be convenient to use p$^\dagger$q to generate the corresponding NumPy einsum code for us. The following code does that.

from pdaggerq.parser import contracted_strings_to_tensor_terms

print('')
print('# aaaa terms')
print('')

spin_labels = {
    'a' : 'a',
    'b' : 'a',
    'i' : 'a',
    'j' : 'a'
}

terms = pq.strings(spin_labels = spin_labels)
terms = contracted_strings_to_tensor_terms(terms)
for term in terms:
    print("#\t", term)
    print("%s" % (term.einsum_string(update_val='sigma2_aaaa', output_variables=('a', 'b', 'i', 'j'))))
    print()

print('')
print('# bbbb terms')
print('')

spin_labels = {
    'a' : 'b',
    'b' : 'b',
    'i' : 'b',
    'j' : 'b'
}

terms = pq.strings(spin_labels = spin_labels)
terms = contracted_strings_to_tensor_terms(terms)
for term in terms:
    print("#\t", term)
    print("%s" % (term.einsum_string(update_val='sigma2_bbbb', output_variables=('a', 'b', 'i', 'j'))))
    print()

print('')
print('# abab terms')
print('')

spin_labels = {
    'a' : 'a',
    'b' : 'b',
    'i' : 'a',
    'j' : 'b'
}

terms = pq.strings(spin_labels = spin_labels)
terms = contracted_strings_to_tensor_terms(terms)
for term in terms:
    print("#\t", term)
    print("%s" % (term.einsum_string(update_val='sigma2_abab', output_variables=('a', 'b', 'i', 'j'))))
    print()

# aaaa terms

#	 -0.50 <l,k||l,k>_aaaa*r2_aaaa(a,b,i,j)
sigma2_aaaa += -0.50 * einsum('lklk,abij->abij', g_aaaa[oa, oa, oa, oa], r2_aaaa)

#	 -0.50 <l,k||l,k>_abab*r2_aaaa(a,b,i,j)
sigma2_aaaa += -0.50 * einsum('lklk,abij->abij', g_abab[oa, ob, oa, ob], r2_aaaa)

#	 -0.50 <k,l||k,l>_abab*r2_aaaa(a,b,i,j)
sigma2_aaaa += -0.50 * einsum('klkl,abij->abij', g_abab[oa, ob, oa, ob], r2_aaaa)

#	 -0.50 <l,k||l,k>_bbbb*r2_aaaa(a,b,i,j)
sigma2_aaaa += -0.50 * einsum('lklk,abij->abij', g_bbbb[ob, ob, ob, ob], r2_aaaa)

#	  0.50 <l,k||i,j>_aaaa*r2_aaaa(a,b,l,k)
sigma2_aaaa +=  0.50 * einsum('lkij,ablk->abij', g_aaaa[oa, oa, oa, oa], r2_aaaa)

#	  1.00 P(i,j)*P(a,b)<k,a||c,j>_aaaa*r2_aaaa(c,b,i,k)
contracted_intermediate =  1.00 * einsum('kacj,cbik->abij', g_aaaa[oa, va, va, oa], r2_aaaa)
sigma2_aaaa +=  1.00000 * contracted_intermediate + -1.00000 * einsum('abij->abji', contracted_intermediate)  + -1.00000 * einsum('abij->baij', contracted_intermediate)  +  1.00000 * einsum('abij->baji', contracted_intermediate) 

#	 -1.00 P(i,j)*P(a,b)<a,k||j,c>_abab*r2_abab(b,c,i,k)
contracted_intermediate = -1.00 * einsum('akjc,bcik->abij', g_abab[va, ob, oa, vb], r2_abab)
sigma2_aaaa +=  1.00000 * contracted_intermediate + -1.00000 * einsum('abij->abji', contracted_intermediate)  + -1.00000 * einsum('abij->baij', contracted_intermediate)  +  1.00000 * einsum('abij->baji', contracted_intermediate) 

#	  0.50 <a,b||c,d>_aaaa*r2_aaaa(c,d,i,j)
sigma2_aaaa +=  0.50 * einsum('abcd,cdij->abij', g_aaaa[va, va, va, va], r2_aaaa)


# bbbb terms

#	 -0.50 <l,k||l,k>_aaaa*r2_bbbb(a,b,i,j)
sigma2_bbbb += -0.50 * einsum('lklk,abij->abij', g_aaaa[oa, oa, oa, oa], r2_bbbb)

#	 -0.50 <l,k||l,k>_abab*r2_bbbb(a,b,i,j)
sigma2_bbbb += -0.50 * einsum('lklk,abij->abij', g_abab[oa, ob, oa, ob], r2_bbbb)

#	 -0.50 <k,l||k,l>_abab*r2_bbbb(a,b,i,j)
sigma2_bbbb += -0.50 * einsum('klkl,abij->abij', g_abab[oa, ob, oa, ob], r2_bbbb)

#	 -0.50 <l,k||l,k>_bbbb*r2_bbbb(a,b,i,j)
sigma2_bbbb += -0.50 * einsum('lklk,abij->abij', g_bbbb[ob, ob, ob, ob], r2_bbbb)

#	  0.50 <l,k||i,j>_bbbb*r2_bbbb(a,b,l,k)
sigma2_bbbb +=  0.50 * einsum('lkij,ablk->abij', g_bbbb[ob, ob, ob, ob], r2_bbbb)

#	 -1.00 P(i,j)*P(a,b)<k,a||c,j>_abab*r2_abab(c,b,k,i)
contracted_intermediate = -1.00 * einsum('kacj,cbki->abij', g_abab[oa, vb, va, ob], r2_abab)
sigma2_bbbb +=  1.00000 * contracted_intermediate + -1.00000 * einsum('abij->abji', contracted_intermediate)  + -1.00000 * einsum('abij->baij', contracted_intermediate)  +  1.00000 * einsum('abij->baji', contracted_intermediate) 

#	  1.00 P(i,j)*P(a,b)<k,a||c,j>_bbbb*r2_bbbb(c,b,i,k)
contracted_intermediate =  1.00 * einsum('kacj,cbik->abij', g_bbbb[ob, vb, vb, ob], r2_bbbb)
sigma2_bbbb +=  1.00000 * contracted_intermediate + -1.00000 * einsum('abij->abji', contracted_intermediate)  + -1.00000 * einsum('abij->baij', contracted_intermediate)  +  1.00000 * einsum('abij->baji', contracted_intermediate) 

#	  0.50 <a,b||c,d>_bbbb*r2_bbbb(c,d,i,j)
sigma2_bbbb +=  0.50 * einsum('abcd,cdij->abij', g_bbbb[vb, vb, vb, vb], r2_bbbb)


# abab terms

#	 -0.50 <l,k||l,k>_aaaa*r2_abab(a,b,i,j)
sigma2_abab += -0.50 * einsum('lklk,abij->abij', g_aaaa[oa, oa, oa, oa], r2_abab)

#	 -0.50 <l,k||l,k>_abab*r2_abab(a,b,i,j)
sigma2_abab += -0.50 * einsum('lklk,abij->abij', g_abab[oa, ob, oa, ob], r2_abab)

#	 -0.50 <k,l||k,l>_abab*r2_abab(a,b,i,j)
sigma2_abab += -0.50 * einsum('klkl,abij->abij', g_abab[oa, ob, oa, ob], r2_abab)

#	 -0.50 <l,k||l,k>_bbbb*r2_abab(a,b,i,j)
sigma2_abab += -0.50 * einsum('lklk,abij->abij', g_bbbb[ob, ob, ob, ob], r2_abab)

#	  0.50 <l,k||i,j>_abab*r2_abab(a,b,l,k)
sigma2_abab +=  0.50 * einsum('lkij,ablk->abij', g_abab[oa, ob, oa, ob], r2_abab)

#	  0.50 <k,l||i,j>_abab*r2_abab(a,b,k,l)
sigma2_abab +=  0.50 * einsum('klij,abkl->abij', g_abab[oa, ob, oa, ob], r2_abab)

#	 -1.00 <a,k||c,j>_abab*r2_abab(c,b,i,k)
sigma2_abab += -1.00 * einsum('akcj,cbik->abij', g_abab[va, ob, va, ob], r2_abab)

#	 -1.00 <k,b||c,j>_abab*r2_aaaa(c,a,i,k)
sigma2_abab += -1.00 * einsum('kbcj,caik->abij', g_abab[oa, vb, va, ob], r2_aaaa)

#	  1.00 <k,b||c,j>_bbbb*r2_abab(a,c,i,k)
sigma2_abab +=  1.00 * einsum('kbcj,acik->abij', g_bbbb[ob, vb, vb, ob], r2_abab)

#	  1.00 <k,a||c,i>_aaaa*r2_abab(c,b,k,j)
sigma2_abab +=  1.00 * einsum('kaci,cbkj->abij', g_aaaa[oa, va, va, oa], r2_abab)

#	 -1.00 <a,k||i,c>_abab*r2_bbbb(c,b,j,k)
sigma2_abab += -1.00 * einsum('akic,cbjk->abij', g_abab[va, ob, oa, vb], r2_bbbb)

#	 -1.00 <k,b||i,c>_abab*r2_abab(a,c,k,j)
sigma2_abab += -1.00 * einsum('kbic,ackj->abij', g_abab[oa, vb, oa, vb], r2_abab)

#	  0.50 <a,b||c,d>_abab*r2_abab(c,d,i,j)
sigma2_abab +=  0.50 * einsum('abcd,cdij->abij', g_abab[va, vb, va, vb], r2_abab)

#	  0.50 <a,b||d,c>_abab*r2_abab(d,c,i,j)
sigma2_abab +=  0.50 * einsum('abdc,dcij->abij', g_abab[va, vb, va, vb], r2_abab)

Here, 'oa', 'ob', 'va', and 'vb' are slices over occupied and virtual orbital labels. Below, we implement these expressions and evaluates the third-order contribution to the energy. First, we need to evaluate the first-order amplitudes.

# construct first-order amplitudes

# aaaa
r2_aaaa = np.zeros([nmo-nalpha, nmo-nalpha, nalpha, nalpha])
for i in range (0, nalpha):
    for j in range (0, nalpha):
        for a in range (nalpha, nmo):
            for b in range (nalpha, nmo):
                denom = epsilon_a[i] + epsilon_a[j] - epsilon_a[a] - epsilon_a[b]
                r2_aaaa[a-nalpha, b-nalpha, i, j] = g_aaaa[i, j, a, b] / denom

# bbbb
r2_bbbb = np.zeros([nmo-nbeta, nmo-nbeta, nbeta, nbeta])
for i in range (0, nbeta):
    for j in range (0, nbeta):
        for a in range (nbeta, nmo):
            for b in range (nbeta, nmo):
                denom = epsilon_b[i] + epsilon_b[j] - epsilon_b[a] - epsilon_b[b]
                r2_bbbb[a-nbeta, b-nbeta, i, j] = g_bbbb[i, j, a, b] / denom


# abab
r2_abab = np.zeros([nmo-nalpha, nmo-nbeta, nalpha, nbeta])
for i in range (0, nalpha):
    for j in range (0, nbeta):
        for a in range (nalpha, nmo):
            for b in range (nbeta, nmo):
                denom = epsilon_a[i] + epsilon_b[j] - epsilon_a[a] - epsilon_b[b]
                r2_abab[a-nalpha, b-nbeta, i, j] = g_abab[i, j, a, b] / denom

Now, we can evaluate

$$\begin{align} \sigma_{ab,ij} = \langle \Phi_{ij}^{ab} | \hat{v} |\psi_0^{(1)} \rangle = \frac{1}{4} \sum_{kl,cd} \langle \Phi_{ij}^{ab} | \hat{v} | \Phi_{kl}^{cd} \rangle c^{(1)}_{cd,kl} \end{align}$$

for each spin case using the code generated by p$^\dagger$q. For that code to function correctly, we must add a few lines to define the slices 'oa', etc. and to build each of the containers, sigma2_aaaa, etc.

from numpy import einsum

# define slices

oa = slice(None, nalpha)
ob = slice(None, nbeta)
va = slice(nalpha, None)
vb = slice(nbeta, None)

# evaluate < ab,ij | v | psi(1) >

# aaaa terms
sigma2_aaaa = np.zeros_like(r2_aaaa)

#	 -0.50 <l,k||l,k>_aaaa*r2_aaaa(a,b,i,j)
sigma2_aaaa += -0.50 * einsum('lklk,abij->abij', g_aaaa[oa, oa, oa, oa], r2_aaaa)

#	 -0.50 <l,k||l,k>_abab*r2_aaaa(a,b,i,j)
sigma2_aaaa += -0.50 * einsum('lklk,abij->abij', g_abab[oa, ob, oa, ob], r2_aaaa)

#	 -0.50 <k,l||k,l>_abab*r2_aaaa(a,b,i,j)
sigma2_aaaa += -0.50 * einsum('klkl,abij->abij', g_abab[oa, ob, oa, ob], r2_aaaa)

#	 -0.50 <l,k||l,k>_bbbb*r2_aaaa(a,b,i,j)
sigma2_aaaa += -0.50 * einsum('lklk,abij->abij', g_bbbb[ob, ob, ob, ob], r2_aaaa)

#	  0.50 <l,k||i,j>_aaaa*r2_aaaa(a,b,l,k)
sigma2_aaaa +=  0.50 * einsum('lkij,ablk->abij', g_aaaa[oa, oa, oa, oa], r2_aaaa)

#	  1.00 P(i,j)*P(a,b)<k,a||c,j>_aaaa*r2_aaaa(c,b,i,k)
contracted_intermediate =  1.00 * einsum('kacj,cbik->abij', g_aaaa[oa, va, va, oa], r2_aaaa)
sigma2_aaaa +=  1.00000 * contracted_intermediate + -1.00000 * einsum('abij->abji', contracted_intermediate)  + -1.00000 * einsum('abij->baij', contracted_intermediate)  +  1.00000 * einsum('abij->baji', contracted_intermediate) 

#	 -1.00 P(i,j)*P(a,b)<a,k||j,c>_abab*r2_abab(b,c,i,k)
contracted_intermediate = -1.00 * einsum('akjc,bcik->abij', g_abab[va, ob, oa, vb], r2_abab)
sigma2_aaaa +=  1.00000 * contracted_intermediate + -1.00000 * einsum('abij->abji', contracted_intermediate)  + -1.00000 * einsum('abij->baij', contracted_intermediate)  +  1.00000 * einsum('abij->baji', contracted_intermediate) 

#	  0.50 <a,b||c,d>_aaaa*r2_aaaa(c,d,i,j)
sigma2_aaaa +=  0.50 * einsum('abcd,cdij->abij', g_aaaa[va, va, va, va], r2_aaaa)


# bbbb terms
sigma2_bbbb = np.zeros_like(r2_bbbb)

#	 -0.50 <l,k||l,k>_aaaa*r2_bbbb(a,b,i,j)
sigma2_bbbb += -0.50 * einsum('lklk,abij->abij', g_aaaa[oa, oa, oa, oa], r2_bbbb)

#	 -0.50 <l,k||l,k>_abab*r2_bbbb(a,b,i,j)
sigma2_bbbb += -0.50 * einsum('lklk,abij->abij', g_abab[oa, ob, oa, ob], r2_bbbb)

#	 -0.50 <k,l||k,l>_abab*r2_bbbb(a,b,i,j)
sigma2_bbbb += -0.50 * einsum('klkl,abij->abij', g_abab[oa, ob, oa, ob], r2_bbbb)

#	 -0.50 <l,k||l,k>_bbbb*r2_bbbb(a,b,i,j)
sigma2_bbbb += -0.50 * einsum('lklk,abij->abij', g_bbbb[ob, ob, ob, ob], r2_bbbb)

#	  0.50 <l,k||i,j>_bbbb*r2_bbbb(a,b,l,k)
sigma2_bbbb +=  0.50 * einsum('lkij,ablk->abij', g_bbbb[ob, ob, ob, ob], r2_bbbb)

#	 -1.00 P(i,j)*P(a,b)<k,a||c,j>_abab*r2_abab(c,b,k,i)
contracted_intermediate = -1.00 * einsum('kacj,cbki->abij', g_abab[oa, vb, va, ob], r2_abab)
sigma2_bbbb +=  1.00000 * contracted_intermediate + -1.00000 * einsum('abij->abji', contracted_intermediate)  + -1.00000 * einsum('abij->baij', contracted_intermediate)  +  1.00000 * einsum('abij->baji', contracted_intermediate) 

#	  1.00 P(i,j)*P(a,b)<k,a||c,j>_bbbb*r2_bbbb(c,b,i,k)
contracted_intermediate =  1.00 * einsum('kacj,cbik->abij', g_bbbb[ob, vb, vb, ob], r2_bbbb)
sigma2_bbbb +=  1.00000 * contracted_intermediate + -1.00000 * einsum('abij->abji', contracted_intermediate)  + -1.00000 * einsum('abij->baij', contracted_intermediate)  +  1.00000 * einsum('abij->baji', contracted_intermediate) 

#	  0.50 <a,b||c,d>_bbbb*r2_bbbb(c,d,i,j)
sigma2_bbbb +=  0.50 * einsum('abcd,cdij->abij', g_bbbb[vb, vb, vb, vb], r2_bbbb)


# abab terms
sigma2_abab = np.zeros_like(r2_abab)

#	 -0.50 <l,k||l,k>_aaaa*r2_abab(a,b,i,j)
sigma2_abab += -0.50 * einsum('lklk,abij->abij', g_aaaa[oa, oa, oa, oa], r2_abab)

#	 -0.50 <l,k||l,k>_abab*r2_abab(a,b,i,j)
sigma2_abab += -0.50 * einsum('lklk,abij->abij', g_abab[oa, ob, oa, ob], r2_abab)

#	 -0.50 <k,l||k,l>_abab*r2_abab(a,b,i,j)
sigma2_abab += -0.50 * einsum('klkl,abij->abij', g_abab[oa, ob, oa, ob], r2_abab)

#	 -0.50 <l,k||l,k>_bbbb*r2_abab(a,b,i,j)
sigma2_abab += -0.50 * einsum('lklk,abij->abij', g_bbbb[ob, ob, ob, ob], r2_abab)

#	  0.50 <l,k||i,j>_abab*r2_abab(a,b,l,k)
sigma2_abab +=  0.50 * einsum('lkij,ablk->abij', g_abab[oa, ob, oa, ob], r2_abab)

#	  0.50 <k,l||i,j>_abab*r2_abab(a,b,k,l)
sigma2_abab +=  0.50 * einsum('klij,abkl->abij', g_abab[oa, ob, oa, ob], r2_abab)

#	 -1.00 <a,k||c,j>_abab*r2_abab(c,b,i,k)
sigma2_abab += -1.00 * einsum('akcj,cbik->abij', g_abab[va, ob, va, ob], r2_abab)

#	 -1.00 <k,b||c,j>_abab*r2_aaaa(c,a,i,k)
sigma2_abab += -1.00 * einsum('kbcj,caik->abij', g_abab[oa, vb, va, ob], r2_aaaa)

#	  1.00 <k,b||c,j>_bbbb*r2_abab(a,c,i,k)
sigma2_abab +=  1.00 * einsum('kbcj,acik->abij', g_bbbb[ob, vb, vb, ob], r2_abab)

#	  1.00 <k,a||c,i>_aaaa*r2_abab(c,b,k,j)
sigma2_abab +=  1.00 * einsum('kaci,cbkj->abij', g_aaaa[oa, va, va, oa], r2_abab)

#	 -1.00 <a,k||i,c>_abab*r2_bbbb(c,b,j,k)
sigma2_abab += -1.00 * einsum('akic,cbjk->abij', g_abab[va, ob, oa, vb], r2_bbbb)

#	 -1.00 <k,b||i,c>_abab*r2_abab(a,c,k,j)
sigma2_abab += -1.00 * einsum('kbic,ackj->abij', g_abab[oa, vb, oa, vb], r2_abab)

#	  0.50 <a,b||c,d>_abab*r2_abab(c,d,i,j)
sigma2_abab +=  0.50 * einsum('abcd,cdij->abij', g_abab[va, vb, va, vb], r2_abab)

#	  0.50 <a,b||d,c>_abab*r2_abab(d,c,i,j)
sigma2_abab +=  0.50 * einsum('abdc,dcij->abij', g_abab[va, vb, va, vb], r2_abab)

Now, the last ingredient we need in order to construct the second-order amplitudes is the first-order energy. We derived spin-orbital expressions for $E_0^{(1)}$ above, but since we are working in a spatial orbital basis, we can let p$^\dagger$q give us code that includes spin labels.

pq = pdaggerq.pq_helper('fermi')

pq.add_operator_product(1.0, ['v'])
pq.simplify()

terms = pq.strings(spin_labels = {})
terms = contracted_strings_to_tensor_terms(terms)
for term in terms:
    print("#\t", term)
    print("%s" % (term.einsum_string(update_val='e1')))
    print()

#	 -0.50 <j,i||j,i>_aaaa
e1 += -0.50 * einsum('jiji', g_aaaa[oa, oa, oa, oa])

#	 -0.50 <j,i||j,i>_abab
e1 += -0.50 * einsum('jiji', g_abab[oa, ob, oa, ob])

#	 -0.50 <i,j||i,j>_abab
e1 += -0.50 * einsum('ijij', g_abab[oa, ob, oa, ob])

#	 -0.50 <j,i||j,i>_bbbb
e1 += -0.50 * einsum('jiji', g_bbbb[ob, ob, ob, ob])

The following code implements these expressions to evaluate the first-order correction to the energy.

# first-order energy, E(1) = < 0 | v | 0 >

e1 = 0.0 

#	 -0.50 <j,i||j,i>_aaaa
e1 += -0.50 * einsum('jiji', g_aaaa[oa, oa, oa, oa])

#	 -0.50 <j,i||j,i>_abab
e1 += -0.50 * einsum('jiji', g_abab[oa, ob, oa, ob])

#	 -0.50 <i,j||i,j>_abab
e1 += -0.50 * einsum('ijij', g_abab[oa, ob, oa, ob])

#	 -0.50 <j,i||j,i>_bbbb
e1 += -0.50 * einsum('jiji', g_bbbb[ob, ob, ob, ob])

Now, we are ready to construct the second-order amplitudes and evaluate $E^{(3)}_0$.

# construct second-order amplitudes

# aaaa
amp2_aaaa = sigma2_aaaa - e1 * r2_aaaa
for i in range (0, nalpha):
    for j in range (0, nalpha):
        for a in range (nalpha, nmo):
            for b in range (nalpha, nmo):
                denom = epsilon_a[i] + epsilon_a[j] - epsilon_a[a] - epsilon_a[b]
                amp2_aaaa[a-nalpha, b-nalpha, i, j] /= denom

# bbbb
amp2_bbbb = sigma2_bbbb - e1 * r2_bbbb
for i in range (0, nbeta):
    for j in range (0, nbeta):
        for a in range (nbeta, nmo):
            for b in range (nbeta, nmo):
                denom = epsilon_b[i] + epsilon_b[j] - epsilon_b[a] - epsilon_b[b]
                amp2_bbbb[a-nbeta, b-nbeta, i, j] /= denom


# abab
amp2_abab = sigma2_abab - e1 * r2_abab
for i in range (0, nalpha):
    for j in range (0, nbeta):
        for a in range (nalpha, nmo):
            for b in range (nbeta, nmo):
                denom = epsilon_a[i] + epsilon_b[j] - epsilon_a[a] - epsilon_b[b]
                amp2_abab[a-nalpha, b-nbeta, i, j] /= denom

# third order energy:

e3 = 0.0
e3 += 0.25 * np.einsum('abij,abij->', g_aaaa[va, va, oa, oa], amp2_aaaa)
e3 += 0.25 * np.einsum('abij,abij->', g_bbbb[vb, vb, ob, ob], amp2_bbbb)
e3 += np.einsum('abij,abij->', g_abab[va, vb, oa, ob], amp2_abab)

e_mp3 = psi4.energy('mp3')
e3_psi4 = e_mp3 - e_mp2

print('')
print('my E(3) =   %20.12f Eh' % (e3))
print('psi4 E(3) = %20.12f Eh' % (e3_psi4))
print('')

my E(3) =        -0.007965387470 Eh
psi4 E(3) =      -0.007965387470 Eh

We can see that our third-order energy is numerically equivalent to that from Psi4. What a day!

© 2026. All rights reserved.