Perturbation Theory¶

In the last notebook, we used the variation theorem to find approximate solutions to the Schrödinger equation. An alterative approach is based on perturbation theory. This notebook outlines the fundamental principles of time-independent and time-dependent perturbation theories.

Time-Independent Perturbation Theory¶

Time-independent perturbation theory is an approach for finding approximate solutions to the time-independent Schrödinger equation that complements the variational approach discussed in the last notebook. The basic idea is that the Hamiltonian for a system can be partitioned into a "simple" Hamiltonian for which solutions can be determined and a "small" perturbation to that simple Hamiltonian. We have

$$\begin{align} \hat{H} = \hat{H}_0 + \lambda \hat{H}^\prime \end{align}$$

Here, $H_0$ is the zeroth-order Hamiltonian, which is the simple part of the problem and $\hat{H}^\prime$ is the operator representing the perturbation. We have also introduced a parameter, $\lambda,$ which can have any value on the interval $[0, 1].$ In the limit that $\lambda \to 0,$ we recover the simple Hamiltonian. In the limit that $\lambda \to 1,$ we recover the full Hamiltonian. We assume that we can determine the eigenfunctions of $\hat{H}_0$ that satisfy

$$\begin{align} \hat{H}_0 | \psi_n^{(0)}\rangle = E_n^{(0)}| \psi_n^{(0)}\rangle \end{align}$$

where $E_n^{(0)}$ and $\psi_n^{(0)}$ are called the zeroth-order energies and wave functions, respectively.

Non-Degenerate Perturbation Theory¶

The next steps differ depending on whether or not any of the eigenvalues of $\hat{H}_0$ are degenerate. We will consider only the case that the eigenvalues are all non-degenerate; the resulting perturbation theory is thus called non-denerate perturbation theory.

Consider a wave function, $|\psi_n\rangle$, that satisfies the time-independent Schrödinger equation,

$$\begin{align} \hat{H} |\psi_n\rangle = E_n |\psi_n\rangle \end{align}$$

$$\begin{align} (\hat{H}_0 + \lambda \hat{H}^\prime) |\psi_n\rangle = E_n |\psi_n\rangle \end{align}$$

where the subscript $n$ refers to a particular state. We now expand $|\psi_n\rangle$ as a Taylor series in the perturbation strength, $\lambda$:

$$\begin{align} |\psi_n\rangle = \left . |\psi_n\rangle \right |_{\lambda = 0} + \lambda \left . \frac{\partial |\psi_n\rangle}{\partial \lambda} \right |_{\lambda = 0} + \lambda ^2\frac{1}{2!}\left . \frac{\partial^2 |\psi_n\rangle}{\partial \lambda^2} \right |_{\lambda = 0} + \lambda ^3 \frac{1}{3!} \left . \frac{\partial^3 |\psi_n\rangle}{\partial \lambda^3} \right |_{\lambda = 0} + ... \end{align}$$

We can simply our notation by first noting that

$$\begin{align} \left . |\psi_n\rangle \right |_{\lambda = 0} = |\psi_n^{(0)}\rangle \end{align}$$

Next, we define the $p$-th order correction to the wave function in terms of the $p$-th order derivative of the wave function, evaluated at $\lambda = 0$, i.e.,

$$\begin{align} \frac{1}{p!} \left . \frac{\partial^p |\psi_n\rangle}{\partial \lambda^p} \right |_{\lambda = 0} = |\psi_n^{(p)}\rangle \end{align}$$

Now, the full wave function can be expressed 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}$$

The energy for state $n$ is expanded in a similar way:

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

where $E_n^{(p)}$ is referred to as the $p$-th order correction to the energy. In what follows, we will derive expressions for the $p$-th order corrections to the wave function and energy, which are typically determined order-by-order, starting with $p = 1$.

If the corrections to the wave function are determined order-by-order, how is the total wave function normalized? Here, we use a special type of normalization condition, called intermediate normalization in which the overlap of the full wave function with the zeroth-order wave function is one:

$$\begin{align} \langle \psi_n^{(0)} | \psi_n \rangle = 1 \end{align}$$

If we expand $|\psi_n\rangle$ in powers of $\lambda$, we have

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

If we choose $|\psi_n^{(0)}\rangle$ to be normalized, then we have

$$\begin{align} \langle \psi_n^{(0)} | \psi_n^{(0)}\rangle = 1 \end{align}$$

and

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

This expression must be satisfied for all powers of $\lambda$. As a result, we can conclude that

$$\begin{align} \langle \psi_n^{(0)} | \psi_n^{(p)}\rangle = 0\text{, }\forall\text{ }p > 0 \end{align}$$

In other words, the zeroth-order wave function for state $n$, $|\psi_n^{(0)}\rangle$, is orthogonal to all corrections to the wave function, $|\psi_n^{(p)}\rangle$. Note that the reason we can assume that $|\psi_n^{(0)}\rangle$ is normalized is that it is an eigenfunction of a Hermitian operator ($\hat{H}_0$) and is therefore normalizable. For the remainder of this discussion, we will assume that the zeroth-order wave functions form an orthonormal set, i.e.,

$$\begin{align} \langle \psi_n^{(0)}|\psi_m^{(0)}\rangle = \delta_{mn} \end{align}$$

Let us go back to the Schrödinger equation, with $|\psi_n\rangle$ and $E_n$ expanded in powers of $\lambda$:

$$\begin{align} (\hat{H}_0 + \lambda \hat{H}^\prime) ( |\psi_0\rangle + \lambda |\psi_n^{(1)}\rangle + ...) = ( E_n^{(0)} + \lambda E_n^{(1)} + ...)( |\psi_0\rangle + \lambda |\psi_n^{(1)}\rangle + ...) \end{align}$$

The right- and left-hand sides of this equation must be equal for all possible values of $\lambda$, which implies that terms on either side having the same power of $\lambda$ must be equivalent. Each order in $\lambda$ defines a different equation involving only a few of the corrections to the wave function and energy. For up to $\lambda = 2$, we have

$$\begin{align} \lambda^0:~~& \hat{H}_0 |\psi_n^{(0)}\rangle = E_n^{(0)} |\psi_n^{(0)}\rangle \\ \lambda^1:~~& \hat{H}_0 |\psi_n^{(1)}\rangle + \hat{H}^\prime | \psi_n^{(0)} \rangle = E_n^{(0)} |\psi_n^{(1)}\rangle + E_n^{(1)} |\psi_n^{(0)}\rangle \\ \lambda^2:~~& \hat{H}_0 |\psi_n^{(2)}\rangle + \hat{H}^\prime | \psi_n^{(1)} \rangle = E_n^{(0)} |\psi_n^{(2)}\rangle + E_n^{(1)} |\psi_n^{(1)}\rangle + E_n^{(2)} |\psi_n^{(0)}\rangle\\ ... \end{align}$$

The first of these equations is simply the zeroth-order Schrödinger equation.

The First-Order Correction to the Energy¶

What is the first-order correction to the energy for state $n$, $E_n^{(1)}$? We can determine this quantity using the first-order Schrödinger equation,

$$\begin{align} \hat{H}_0 |\psi_n^{(1)}\rangle + \hat{H}^\prime | \psi_n^{(0)} \rangle = E_n^{(0)} |\psi_n^{(1)}\rangle + E_n^{(1)} |\psi_n^{(0)}\rangle \end{align}$$

To isolate $E_n^{(1)}$, we can left-multiply both sides of this equation by $\langle \psi_n^{(0)}|$ and integrating each term over all space

$$\begin{align} \langle \psi_n^{(0)}| \hat{H}_0 |\psi_n^{(1)}\rangle + \langle \psi_n^{(0)}| \hat{H}^\prime | \psi_n^{(0)} \rangle = E_n^{(0)} \langle \psi_n^{(0)}|\psi_n^{(1)}\rangle + E_n^{(1)} \langle \psi_n^{(0)}|\psi_n^{(0)}\rangle \end{align}$$

Now, let us inspect each term and play my favorite game: "zero, one, or something else?!".

$\langle \psi_n^{(0)}| \hat{H}_0 |\psi_n^{(1)}\rangle = 0$. "How can this be?" you might ask. Well, two pieces of information are key. First, $|\psi_n^{(0)}\rangle$ is an eigenfunction of $\hat{H}_0$. Second, $|\psi_n^{(0)}\rangle$ is orthogonal to any correction to the wave function for state $n$ (e.g., $|\psi_n^{(1)}\rangle$). More explicitly, $$\begin{align} \langle \psi_n^{(0)}| \hat{H}_0 |\psi_n^{(1)}\rangle & = \langle \psi_n^{(1)}| \hat{H}_0 |\psi_n^{(0)}\rangle ^* \\ & = E_n^{(0)} \langle \psi_n^{(1)} | \psi_n^{(0)}\rangle^* \\ & = 0 \end{align}$$
$\langle \psi_n^{(0)}| \hat{H}^\prime | \psi_n^{(0)} \rangle = $ something else.
$\langle \psi_n^{(0)}|\psi_n^{(1)}\rangle = 0$ because, as mentioned above, these functions are orthogonal.
$\langle \psi_n^{(0)}|\psi_n^{(0)}\rangle = 1$ because the zeroth-order wave functions are normalized.

Thus, to first-order in perturbation theory,

$$\begin{align} E_n = E_n^{(0)} + \langle \psi_n^{(0)}| \hat{H}^\prime | \psi_n^{(0)} \rangle \end{align}$$

Now, let us evaluate the ground-state energy (to first order) of a one-dimensional anharmonic oscillator of mass $m$ with the following Hamiltonian

$$\begin{align} \hat{H} = -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + \frac{1}{2} k x^2 + c x^3 + d x^4 \end{align}$$

Here, $x$ reprents the displacement of the oscillator from its equilibrium position, and $k,$ $c,$ and $d$ are constants. We can partition this Hamiltonian into a zeroth-order part that looks like the Hamiltonian for the quantum harmonic oscillator (QHO) problem,

$$\begin{align} \hat{H}_0 = -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + \frac{1}{2} k x^2 \end{align}$$

and a perturbation

$$\begin{align} \hat{H}^\prime = c x^3 + d x^4 \end{align}$$

The zeroth-order wave functions have the form

$$\begin{align} |\psi_n^{(0)}\rangle = N_n e^{-\alpha x^2/2} H_n(\alpha^{1/2}x) \end{align}$$

where $\alpha = \frac{\omega m}{\hbar}$, $\omega = 2\pi \nu$, $\nu$ is the frequency of the oscillations, $N_n$ is a normalization constant, and $H_n$ is a Hermite polynomial. The zeroth-order energies are

$$\begin{align} E_n^{(0)} = \left (n + \frac{1}{2} \right ) \hbar \omega \end{align}$$

To first order, the ground-state energy is

$$\begin{align} E_0 = E_0^{(0)} + \langle \psi_0^{(0)} | \hat{H}^\prime | \psi_0^{(0)}\rangle \end{align}$$

Given that $E_0^{(0)} = \frac{1}{2}\hbar \omega$ and $|\psi_0^{(0)}\rangle = \left ( \frac{\alpha}{\pi}\right ) ^{1/4}e^{-\alpha x^2/2} $, we have

$$\begin{align} E_0 &= \frac{1}{2}\hbar \omega + \int_{-\infty}^{\infty} ( cx^3 + dx^4) \left ( \frac{\alpha}{\pi}\right ) ^{1/2}e^{-\alpha x^2} \\ &= \frac{1}{2}\hbar \omega + d \left ( \frac{\alpha}{\pi}\right ) ^{1/2} \int_{-\infty}^{\infty} x^4 e^{-\alpha x^2} dx \\ &= \frac{1}{2}\hbar \omega + \frac{3d}{4\alpha^2} \end{align}$$

Note that the parameter, $c$, does not appear at first order, but it will appear at higher orders in perturbation theory.

The First-Order Correction to the Wave Function¶

We isolated the first-order correction to the energy by left-multiplying the first-order Schrödinger equation by $\langle \psi_n^{(0)}|$ and integrating over all space. Let us see what happens if instead we project the first-order Schrödinger equation,

$$\begin{align} \hat{H}_0 |\psi_n^{(1)}\rangle + \hat{H}^\prime | \psi_n^{(0)} \rangle = E_n^{(0)} |\psi_n^{(1)}\rangle + E_n^{(1)} |\psi_n^{(0)}\rangle \end{align}$$

onto $\langle \psi_m^{(0)}|$, with $m\neq n$. Doing so gives

$$\begin{align} \langle \psi_m^{(0)}| \hat{H}_0 |\psi_n^{(1)}\rangle + \langle \psi_m^{(0)}| \hat{H}^\prime | \psi_n^{(0)} \rangle = E_n^{(0)} \langle \psi_m^{(0)}|\psi_n^{(1)}\rangle + E_n^{(1)} \langle \psi_m^{(0)}|\psi_n^{(0)}\rangle \end{align}$$

How shall we proceed? You guessed it! It is time for another rousing game of "zero, one, or something else?!".

$\langle \psi_m^{(0)}| \hat{H}_0 |\psi_n^{(1)}\rangle =$ something else. We have $$\begin{align} \langle \psi_m^{(0)}| \hat{H}_0 |\psi_n^{(1)}\rangle & = \langle \psi_n^{(1)}| \hat{H}_0 |\psi_m^{(0)}\rangle ^* \\ & = E_m^{(0)} \langle \psi_n^{(1)} | \psi_m^{(0)}\rangle^* \\ & = E_m^{(0)} \langle \psi_m^{(0)} | \psi_n^{(1)}\rangle \end{align}$$ Note that, while the first-order correction to the wave function for state $n$ will be orthogonal to the zeroth-order wave function for state $n$, we have no such guarantee for the zeroth-order wave function for a different state (in this case, state $m$).
$\langle \psi_m^{(0)}| \hat{H}^\prime | \psi_n^{(0)} \rangle = $ something else.
$\langle \psi_m^{(0)}|\psi_n^{(1)}\rangle = $ something else. As stated above, we have no guarantee that a correction to the wave function for state $n$ is orthogonal to the zeroth-order wave function for a different state, $m$.
$\langle \psi_m^{(0)}|\psi_n^{(0)}\rangle = 0$ because the zeroth-order wave functions form an orthonormal set, and $m \neq n$.

With these considerations, we have

$$\begin{align} \langle \psi_m^{(0)}|\psi_n^{(1)}\rangle = \frac{\langle \psi_m^{(0)}| \hat{H}^\prime | \psi_n^{(0)} \rangle}{E_n^{(0)} - E_m^{(0)}} \end{align}$$

Now, we recall that the zeroth-order wave functions are eigenfunctions of a Hermitian operator, and, they form a complete set. As such, we can expand $|\psi_n^{(1)}\rangle$ as

$$\begin{align} |\psi_n^{(1)}\rangle = \sum_{k\neq n} c_k |\psi_k^{(0)}\rangle \end{align}$$

Note the restriction that $k \neq n$, which is necessary to guaratee that $\langle \psi_n^{(0)}|\psi_n^{(1)}\rangle = 0$. To determine a given expansion coefficient, we must project $|\psi_n^{(1)}\rangle$ onto the corresponding basis function:

$$\begin{align} \langle \psi_m^{(0)} | \psi_n^{(1)} \rangle &= \sum_{k\neq n} c_k\langle \psi_m^{(0)} | \psi_k^{(0)} \rangle \\ &= \sum_{k\neq n} c_k \delta_{mk}\\ &= c_m \end{align}$$

Combined with the result from above, we can see that the expansion coefficients are

$$\begin{align} c_m = \frac{\langle \psi_m^{(0)}| \hat{H}^\prime | \psi_n^{(0)} \rangle}{E_n^{(0)} - E_m^{(0)}} \end{align}$$

so the first-order correction to the wave function is

$$\begin{align} |\psi_n^{(1)}\rangle = \sum_{k\neq n} \frac{\langle \psi_k^{(0)}| \hat{H}^\prime | \psi_n^{(0)} \rangle}{E_n^{(0)} - E_k^{(0)}} | \psi_k^{(0)} \rangle \end{align}$$

From our anharmonic oscillator example, the first-order correction to the wave function for state $n$ would be

$$\begin{align} |\psi_n^{(1)}\rangle = \sum_{k\neq n}^\infty \frac{N_k N_n \int_{-\infty}^{\infty} e^{-\alpha x^2} H_k(\alpha^{1/2}x) H_n(\alpha^{1/2} x) (cx^3 + dx^4)}{\hbar\omega ( n - k )} |\psi_k^{(0)}\rangle \end{align}$$

In this case, we can see that the summation is over an infinite number of states, $k$. Note, however, that this is not generally the case. In Møller-Plesset perturbation theory, for example, evalulating the first-order correction to the wave function does not require an infinite summation.

The Second-Order Correction to the Energy¶

Let's keep this party going. Here, we will derive an expression for the second-order correction to the energy. The second-order energy is particularly important in electronic structure theory. In Møller-Plesset perturbation theory, the second-order correction (MP2) is the first correction to the Hartree-Fock energy.

To obtain the second-order energy, we proceed in a similar manner as above, when we isolated the first-order energy from the first-order Schrödinger equation. The second-order Schrödinger equation is

$$\begin{align} \hat{H}_0 |\psi_n^{(2)}\rangle + \hat{H}^\prime | \psi_n^{(1)} \rangle = E_n^{(0)} |\psi_n^{(2)}\rangle + E_n^{(1)} |\psi_n^{(1)}\rangle + E_n^{(2)} |\psi_n^{(0)}\rangle \end{align}$$

We left-multiply both sides by $\langle \psi_n^{(0)}|$ and integrate over all space, giving

$$\begin{align} \langle \psi_n^{(0)}| \hat{H}_0 |\psi_n^{(2)}\rangle + \langle \psi_n^{(0)} | \hat{H}^\prime | \psi_n^{(1)} \rangle = E_n^{(0)} \langle \psi_n^{(0)}|\psi_n^{(2)}\rangle + E_n^{(1)} \langle \psi_n^{(0)}|\psi_n^{(1)}\rangle + E_n^{(2)} \langle \psi_n^{(0)}|\psi_n^{(0)}\rangle \end{align}$$

Guess what time it is? It is time to play "zero, one, or something else?!". On the left-hand side of this equation, $\langle \psi_n^{(0)}| \hat{H}_0 |\psi_n^{(2)}\rangle = 0$ (you can prove this for yourself). On the right-hand side, only the third term survives because the zeroth order wave function for state $n$ is orthogonal to any correction to the wave function for state $n$. We then have

$$\begin{align} E_n^{(2)} &= \langle \psi_n^{(0)} | \hat{H}^\prime | \psi_n^{(1)}\rangle \end{align}$$

We just showed that

$$\begin{align} |\psi_n^{(1)}\rangle = \sum_{k\neq n} \frac{\langle \psi_k^{(0)}| \hat{H}^\prime | \psi_n^{(0)} \rangle}{E_n^{(0)} - E_k^{(0)}} | \psi_k^{(0)} \rangle \end{align}$$

$$\begin{align} E_n^{(2)} &= \sum_{k\neq n} \langle \psi_n^{(0)} | \hat{H}^\prime | \psi_k^{(0)} \rangle \frac{\langle \psi_k^{(0)}|\hat{H}^\prime | \psi_n^{(0)} \rangle }{E_n^{(0)} - E_k^{(0)}} \\ &= \sum_{k\neq n} \frac{ |\langle \psi_k^{(0)}|\hat{H}^\prime | \psi_n^{(0)} \rangle |^2 }{E_n^{(0)} - E_k^{(0)}} \end{align}$$

Recall the anharmonic oscillator. To second-order the ground-state energy for this system would be

$$\begin{align} E_0 = \frac{1}{2}\hbar \omega + \frac{3d}{4\alpha^2} + \sum_{k\neq 0} \frac{N_k N_0 \left [ \int_{-\infty}^{\infty} e^{-\alpha x^2} H_k(\alpha^{1/2}x) H_0(\alpha^{1/2} x) (cx^3 + dx^4) dx\right ] }{\frac{1}{2}\hbar\omega-\hbar\omega (k + \frac{1}{2})} \end{align}$$

Note that, unlike the first-order energy, the second-order energy for the anharmonic oscillator will depend on both the $c$ and $d$ parameters.

Discussion Points¶

Perturbed energies and wave functions allow us to "mix" character of other zeroth-order states into state $n$. This mixing is necessary because $|\psi_n^{(0)}\rangle$ will not look exactly like $|\psi_n\rangle$.
The $\frac{1}{E_n^{(0)} - E_k^{(0)}}$ term implies that states nearby in energy dominate perturbative corrections for state $n$.
In general, summations can be infinite (e.g., as we saw for the anharmonic oscillator problem). Exceptions include Møller-Plesset perturbation theory.
Higher-order corrections to the energies and wave functions can become extremely complicated. In electronic structure theory, the highest-order corrections one would typically encouter would be fourth or fifth order. As an example, the perturbative triple excitation correction (T) to the coupled-cluster with single and double excitations (CCSD) energy [CCSD(T)], mostly looks like the triples part of the fourth-order energy in Møller-Plesset perturbation theory. It also contains some terms that appear at fifth-order.
If the zeroth-order problem has a continuum of eigenfunctions (e.g., unbound states in the hydrogen atom problem), the these states must be included in the perturbative expansins (in these cases, sums over states become integrals).
The approach outlined above works well for $\hat{H}_0$ with non-degenerate energy levels and when the perturbation is "small." In such cases, $|\psi_n^{(0)}\rangle$ gives a good approximation to $|\psi_n\rangle$.
It is important to recognize that the energy in perturbation theory is not variational, i.e., there is no guarantee that it will be an upper bound to the true ground-state energy.

Møller-Plesset Perturbation theory¶

Møller-Plesset perturbation theory is a form of non-degenerate perturbation theory that is used in electronic structure theory to obtain estimates of the ground-state energy for many-electron systems. The zeroth-order Hamiltonian in Møller-Plesset perturbation theory is the Fock operator, $\hat{F}$, which is an effective one-particle Hamiltonian that arises in Hartree-Fock theory (which is a variational approach and is essentially the mathematical version of molecular orbital theory). Second-order Møller-Plesset perturbation theory (MP2) holds a special place in electronic structure theory because it is the lowest-order of this perturbation theory that improves upon the Hartree-Fock energy. The code below uses the Psi4 package to calculate the energy of molecular hydrogen (H$_2$) as a function of the H-H distance at several levels of theory: Hartree-Fock, full configuration interaction (CI, which is the "exact" variational result within the chosen basis set), and second-, third-, and fourth-order Møller-Plesset perturbation theory (MP2, MP3, and MP4, respectively).

In [33]:

import psi4

import psi4
import numpy as np

# set molecule
mol = psi4.geometry("""
h
h 1 x
symmetry c1
""")   

# set some options for psi4
psi4.set_options({'basis': 'cc-pvdz',
                  'guess_mix': True,
                  'reference': 'rhf'})

# tell psi4 not to print any output to the screen
psi4.core.be_quiet()

# calculate the energy at different H-H separations, x
hf_energy = []
ci_energy = []
mp2_energy = []
mp3_energy = []
mp4_energy = []
dx = 0.1
x = np.arange(0.3, 8.1, dx)
for i in range (0, len(x)):
    mol.x = x[i]
    
    # full ci
    en = psi4.energy('fci')
    ci_energy.append(en)

    # hartree-fock and MPn
    en = psi4.energy('mp4')
    ehf = psi4.variable("HF TOTAL ENERGY")
    hf_energy.append(ehf)
    mp2_energy.append(ehf + psi4.variable("MP2 CORRELATION ENERGY"))
    mp3_energy.append(ehf + psi4.variable("MP3 CORRELATION ENERGY"))
    mp4_energy.append(ehf + psi4.variable("MP4 CORRELATION ENERGY"))

In [34]:

import matplotlib.pyplot as plt

plt.figure()
plt.plot(x, hf_energy, label='Hartree-Fock')
plt.plot(x, ci_energy, label='configuration interaction')
plt.plot(x, mp2_energy, label='MP2')
plt.plot(x, mp3_energy, label='MP3')
plt.plot(x, mp4_energy, label='MP4')

plt.xlim(x[0], x[-1])

plt.ylabel(r'energy (E$_h$)')
plt.xlabel('H-H distance (Å)')
plt.legend()
plt.show()

No description has been provided for this image

The Hartree-Fock approach is variational, so the Hartree-Fock energy is always an upper-bound to the full CI energy. We can see that MP2, MP3, and MP4 improve upon the Hartree-Fock result near equilibrium but begin to give nonsense results at stretched H--H distances. In particular, the MPn energy can be lower than the "exact" variational approach (recall that perturbation theory is not a variational theory). Why is this happening? Like any non-degenerate perturbation theory, Møller-Plesset perturbation theory only works well when the zeroth-order energies are non-degenerate. Recall the expression for the second-order energy,

$$\begin{align} E_n^{(2)} = \sum_{k\neq n} \frac{ |\langle \psi_k^{(0)}|\hat{H}^\prime | \psi_n^{(0)} \rangle |^2 }{E_n^{(0)} - E_k^{(0)}} \end{align}$$

It is clear that this energy could be poorly behaved if $E_n^{(0)} - E_k^{(0)}$ becomes small. In the case of MPn, the zeroth-order energies are sums of orbital energies (the negative of which can be interpreted as the energy required to remove an electron from that orbital). In bond breaking processes, the bonding and antibonding orbital energies approach one another, and we run the risk of encountering small energy denominators in the MP2 energy expression. This issue becomes even worse for higher-order corrections ($E_n^{(3)}$, $E_n^{(4)}$, etc.) because the relevant expressions contain products of these energy differences. Degeneracies and near degeneracies in orbital energies are commonly encountered in chemical applications aside from bond-breaking processes. As such, extreme care should be taken when using perturbation theory in practice.

In [ ]: