Møller-Plesset perturbation theory

Tue 27 August 2019 | tags: blogMP2

Introduction

In the context of theoretical chemistry perturbation theory refers to Rayleigh-Schrödinger perturbation theory. The idea is to divide the Hamilton-operator into two parts, where the solution for one part is known and the second part has a small contribution.

The first part is called the unperturbed part and the second part the perturbation. An example is an molecule in an external electric field. The molecule is than described with Hartree Fock (HF), where the contribution of the external field is described by the perturbation.

\begin{equation} \hat{H}_{0} \psi_{0} = E_{0} \psi_{0} \end{equation}
\begin{equation} E = \langle{\psi_0}| {\hat{H}_{0}}| {\psi_0} \rangle + \langle{\psi_0}|{\hat{V}}|{\psi_0}\rangle \end{equation}

The idea of the perturbation theory is similar to a Taylor expansion of a function. But \({\hat{V}}\) is not a function, but an operator, and so an auxiliary parameter \(\lambda\) has to be introduced.

\begin{equation} {\hat{V}} \rightarrow \lambda {\hat{V}} \end{equation}

It represents the strength of the perturbation and has only two physical meaningful values zero for no perturbation and one for full perturbation. With this adjustments the full Hamiltonian depends on \(\lambda\)

\begin{equation} \hat{H}(\lambda) = \hat{H}_{0} + \lambda \hat{V} \end{equation}

and the energy and wave function can be expressed in a series expansion of \(\lambda\)

\begin{equation} E_n(\lambda) = E_{n}^{(0)} + \sum_k^{\infty} \lambda^k E_n^{(k)} \end{equation}
\begin{equation} \psi_n(\lambda) = \psi_{n}^{(0)} + \sum_k^{\infty} \lambda^k \psi_n^{(k)} \end{equation}

The key idea is that the eigen-values and -functions continuous grow with \(\lambda\), which is only fulfilled for small perturbations.\ If one plugs in the expansion into the Ansatz for the Schrödinger equation one obtains

\begin{equation} \begin{split} & \left( \hat{H}_0 + \lambda \hat{V} \right) \left( \psi_n^{(0)} + \lambda^1 \psi_n^{(1)} + \lambda^2 \psi_n^{(2)} + ... \right) \\ & = \left( E_n^{(0)} + \lambda^1 E_n^{(1)} + \lambda^2 E_n^{(2)} + ... \right) \left( \psi_n^{(0)} + \lambda^1 \psi_n^{(1)} + \lambda^2 \psi_n^{(2)} + ... \right) \end{split} \end{equation}

This equation must be valid for all values of \(\lambda\), which is only the case if it is satisfied for each power of \(\lambda\)

\begin{equation} \label{eq:eigval} \begin{split} & \hat{H}_0 \psi_n^{(0)} = E_n^{(0)} \psi_n^{(0)} \\ & \hat{H}_0 \psi_n^{(1)} + \hat{V} \psi_n^{(0)} = E_n^{(0)} \psi_n^{(1)} + E_n^{(1)} \psi_n^{(0)} \\ & \hat{H}_0 \psi_n^{(2)} + \hat{V} \psi_n^{(1)} = E_n^{(0)} \psi_n^{(2)} + E_n^{(1)} \psi_n^{(1)} + E_n^{(2)} \psi_n^{(0)} \end{split} \end{equation}

All these equations have to be solved successively up to the desired order of perturbation. To make the formalism a bit easier one uses in general the so called intermediate normalization

\begin{equation} \langle{\psi_n^{(0)}}|{\psi_n(\lambda)}\rangle = 1 \end{equation}
\begin{equation} \langle{\psi_n^{(0)}}|{\psi_n^{(k)}}\rangle = 0 \end{equation}

with \(k>0\). In this way higher order contributions to an eigenstate are orthogonal to the zero order contribution.

The energy contribution of first order is obtained after multiplication of the (2nd) eigenvalue equation with \(\langle{\psi_n^{(0)}}|\):

\begin{equation} \langle{\psi_n^{(0)}}|{\hat{V}}|{\psi_n^{(0)}}\rangle = E_n^{(1)} \end{equation}

The wavefunction (WF) of zero order is known and we develop the first order WF in the orthonormal basis of the zero order states \(\psi_m^{(0)}\)

\begin{equation} \psi_n^{(1)} = \sum_{m\ne n} c_{mn} \psi_m^{(0)} \end{equation}

Plugging in this Ansatz for our wave function and multiplying with \(\langle{\psi_j^{(0)}}|\) gives

\begin{equation} \psi_n^{(1)} = \sum_{m\ne n}\frac{\langle{\psi_n^{(0)}}|{\hat{V}}|{\psi_n^{(0)}}\rangle}{E_n^{(0)} - E_m^{(0)}} \psi_m^{(0)} \end{equation}

With the first order WF the second order energy contribution can be calculated as:

\begin{equation} E_n^{(2)} = \langle{\psi_n^{(0)}}|{\hat{V}}|{\psi_n^{(1)}}\rangle = \psi_n^{(1)} = \sum_{m\ne n}\frac{\left|\langle{\psi_n^{(0)}}|{\hat{V}}|{\psi_n^{(0)}}\rangle\right|^{2}} {E_n^{(0)} - E_m^{(0)}} \end{equation}

Møller-Plesset perturbation theory

Perturbation theory is a general tool to add small corrections to a solved system. Møller-Plesset perturbation theory is RS perturbation theory applied to dynamical electron correlation. The reference system is the solved HF problem. So the unperturbed Hamiltonian is defined as the Fock operator

\begin{equation} \hat{H}_0 = \sum_i \hat{F}_i = \sum_i \hat{h}_i + \sum_{ij} \left( \hat{J}_{ij} - \hat{K}_{ij} \right) \end{equation}

and the perturbation is defined as

\begin{equation} \hat{V} = \sum_i \hat{h}_i + \sum_{ij} \frac{1}{2r_{ij}} - \hat{H}_{0} = \sum_{ij} \frac{1}{2r_{ij}} - \left( \hat{J}_{ij} - \hat{K}_{ij} \right) \end{equation}

With this definition it is obvious that the zero order energy is not the HF energy, but

\begin{equation} E_{0}^{(0)} = \langle{\psi_0}|{\hat{H}_0}|{\psi_0}\rangle=\langle{\psi_0}|{\sum_i \hat{F}_i}|{\psi_0}\rangle = \sum_i \varepsilon_i \quad . \end{equation}

The HF energy is obtained with the first order energy correction

\begin{equation} \begin{split} E_0^{(1)} & = \langle{\psi_0}|{\hat{V}}|{\psi_0}\rangle = \frac{1}{2} \sum_{ij} \langle{\psi_0}|{\frac{1}{r_{ij}}}|{\psi_0}\rangle - \sum_{ij} \langle{\psi_0}|{\left(\hat{J}_{ij}-\hat{K}_{ij} \right)}|{\psi_0}\rangle \\ & = -\frac{1}{2}\sum_{ij} \left(\hat{J}_{ij}-\hat{K}_{ij} \right) \quad . \end{split} \end{equation}

The sum of the zero and first order energy contributions equals the HF energy:

\begin{equation} E_0^{(0)} + E_0^{(1)} = \sum_i \varepsilon_i - \frac{1}{2}\sum_{ij} \left(\hat{J}_{ij}-\hat{K}_{ij} \right) = \sum_i \hat{h}_i + \frac{1}{2}\sum_{ij} \left(\hat{J}_{ij}-\hat{K}_{ij} \right) \end{equation}

The energy correction of second order \(E_0^{(2)}\) to the ground state is

\begin{equation} E_0^{(2)} = \sum_{m\ne 0} \frac{\left|\langle{\psi_0^{(0)}}|{\hat{V}}|{\psi_m^{(0)}}\rangle\right|^2}{E_0^{(0)}-E_m^{(0)}} \end{equation}

and contains the interaction of the HF ground state with all excited determinants \(\psi_m^{(0)}\), but only the matrix element with double excited determinants does not vanish. For single excited determinants the expression becomes:

\begin{equation} \langle{\psi_0^{(0)}}|{\hat{V}}|{\psi_1^{(0)}}\rangle = \langle{\psi_0^{(0)}}|{\hat{H}}|{\psi_1^{(0)}}\rangle - \langle{\psi_0^{(0)}}|{\hat{F}}|{\psi_1^{(0)}}\rangle = 0 - \varepsilon_1 \cdot 0 = 0 \end{equation}

The first term vanishes due to the Brillouin theorem and the second term vanishes because of orthogonality. The matrix elements with triple and higher excited determinants vanish, because the perturbation operator is a two electron operator.

If one plugs in the HF wave function one obtains (without derivation) the MP2 energy contribution:

\begin{equation} E_0^{(2)} = \sum_{a<b}^{\rm{vir}}\sum_{i<j}^{\rm{occ}} \frac{\left(\langle{ab}|{ij}\rangle - \langle{ab}|{ji}\rangle \right)^2} {\varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b} \end{equation}

In the literature also the charge cloud notation is often used, where the expression looks like

\begin{equation} E_0^{(2)} = \sum_{a<b}^{\rm{vir}}\sum_{i<j}^{\rm{occ}} \frac{\left( ai|bj \right)^2} {\varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b} \end{equation}

Comments on the computational costs of MP2

The MP2 energy contribution looks quite simple, but is still relatively costly. First the HF problem has to be solved and afterwards the 4 index integrals \(\left( ai|bj \right)\) have to be computed from integrals in the AO basis \(\left( \mu\nu|\lambda\delta \right)\). If one does it in one step the computational costs would scale with the 8\(^{\rm{th}}\) power of the system size \(\mathcal O(N^8)\), because 8 indices occur in the contraction:

\begin{equation} \left( ai|bj \right) = \sum_{\mu\nu\lambda\delta} C_{\mu\nu\lambda\delta}^{aibj} \left( \mu\nu|\lambda\delta \right) \quad . \end{equation}

But with an intelligent summation the costs can be reduced to only \(4\times \mathcal O(N^5)\)

\begin{equation} \left( a\nu|\lambda\delta \right) = \sum_{\mu} C_{\mu}^{a} \left( \mu\nu|\lambda\delta \right) \quad , \end{equation}
\begin{equation} \left( ai|\lambda\delta \right) = \sum_{\nu} C_{\nu}^{i} \left( a\nu|\lambda\delta \right) \quad , \end{equation}
\begin{equation} \left( ai|b\delta \right) = \sum_{\lambda} C_{\lambda}^{b} \left( ai|\lambda\delta \right) \quad , \end{equation}
\begin{equation} \left( ai|bj \right) = \sum_{\delta} C_{\delta}^{j} \left( ai|b\delta \right) \quad . \end{equation}

social